Fitting a simple JAGS model with RSTAN [closed]

Question

I'm trying to fit a simple exponential model for left censored data using RSTAN to replicate something I did in JAGS.

The JAGS model is:

model{
  for(i in 1:N){
    isAboveLOD[i] ~ dinterval(x[i], LOD[i])
    x[i] ~ dexp(lambda)
  }

  lambda ~ dgamma(0.001, 0.001)
}

A full explanation of the simulated data can be found at http://jamescurran.co.nz/2014/06/bayesian-modelling-of-left-censored-data-using-jags/, but the essence of it is that the data had 9691 censored values (out of 10,000) and that the data was originally sampled from an Exp(1.05) distribution, with left censoring performed at log(29).

When I run this model with RSTAN, I essentially get no movement in the chain, and hence the estimate of lambda is way off (0.3 in my first run). I'm sure the error is elementary, but I'd like your advice. My RSTAN code (including data simulation) is given below.

expCode = '
  data {
    int<lower=0> nObs;
    int<lower=0> nCens;
    real<lower=0> yObs[nObs];
    real<lower=0> U;
  }
  parameters {
    real<lower=0, upper=U> yCens[nCens];
    real<lower=0> lambda;
  }
  model {
    for (n in 1:nObs)
      yObs[n] ~ exponential(lambda);
    for (n in 1:nCens)
      yCens[n] ~ exponential(lambda);
    lambda ~ gamma(0.001, 0.001);
  }  
'


set.seed(35202)
x = rexp(10000, rate = 1.05)

## set all the censored values to NA's
x[x < log(29)] = NA

yObs = x[!is.na(x)]
yCens = x[is.na(x)]

dataList = list(nObs = length(yObs), nCens = length(yCens), U = log(29), yObs = yObs)
library(rstan)
fit = stan(model_code = expCode, data = dataList, iter = 1000, chains = 1)

Answer 1

In BUGS/JAGS, the order in which statements are written does not matter. In Stan statements execute in the order in which they are written (see Stan 2.2.0 Reference Manual, pg. 405). Thus your last statement is in the wrong place: lambda is sampled from a gamma distribution, but that happens after the previous statements, so it's sampled 'in the air'.

Furthermore, in Stan

• if no prior is specified for a parameter, it is implicitly given a uniform prior on its support; thus you can define a vague gamma prior for $\lambda$, but that's not needed, real<lower=0> lambda is enough (pg. 8);
• the loop over sampling statements can be vectorized by replacing for (n in 1:N) y[n] ~ ... with the equivalent vectorized form y ~ ... (pg. 9).

If you replace your model block with

  model {
    yObs ~ exponential(lambda);
    yCens ~ exponential(lambda);
  }

you get:

> print(fit, pars="lambda", digits=3)
...
        mean se_mean    sd  2.5%   25%   50%   75% 97.5% n_eff  Rhat
lambda 1.033   0.001 0.015 1.005 1.023 1.034 1.044 1.064   382 0.999
...

EDIT: Please, stop upvoting my answer ;-) There is something wrong in what I've written, but I need a true computer to check what happens here (I'm using a small netbook now). I'll edit my answer as soon as possible.

Answer 2

I am not 100% certain this is the right answer; so take it with a grain of salt. In Stan the orders of statement matters - try putting the prior on lambda before using it. It's possible lambda is just sampling from gamma(0.001, 0.001). Or don't define a prior at all, then you get the uniform distribution (or the improper version of that) on the support as a prior.