censoring model misfit (minimal example)

[Tamas Papp]

Hi,

I ran into a model where the posterior estimate from Stan does not fit
the true value (using simulated data) --- in fact, the region does not
even include it, by a large margin. I have created a minimal example,
posted at https://gist.github.com/tpapp/fadf244cb8c93244412f

The Stan model uses the approach suggested by Bob Carpenter for a more
complicated model
(https://groups.google.com/forum/#!topic/stan-users/PZE0tkwJ2A4).

Before filing an issue, I would apprecite if someone could look at
it. The data generating process is the following:

m <- 1 # boundary
z <- rlnorm(10000) # latent variable
z <- z[z > m] # only kept when above m
N <- length(z)
z_hat <- rlnorm(N, log(z), 0.1) # observed with noise

z is not observed, but z_hat is. The parameter that I am interested in
is m, here = 1, but the stan model (at the link above) estimates it far
below that, moreover, the distribution is not really consistent across
runs (with 1000 iterations).

While it is most likely that I am facing a conceptual issue that I need
to understand, I wonder if I ran into a bug.

Best,

Tamas

[Luc Coffeng]

Hi Tamar,

I think you can only estimate parameter m if you know what the scale of the noise in z_hat is. If you are trying to estimate both the scale of the noise and threshold m, the posterior will entirely depend on what priors you use.

The problem is that the simulated data hold no information on the relative contributions of the noise and the threshold to the marginal variance of the data, unless you have repeated observation of z_hat[ij] | z[j] (which will provide information about the scale of the noise). From your data simulation, it seems you have no repeated observations of z_hat[ij] | z[j], though.

Maybe I'm wrong, but you might want to test it by actually setting the noise term in your model to same value as your data-generating process, or simulate repeated observations z_hat[ij]. Then, threshold m should be identifiable.

Best,

Luc

[Luc Coffeng]

Oops, apologies. I misspelled your name, Tamas!

[Luc Coffeng]

Sorry, I spoke/wrote too soon. Should have read your model code more carefully.

The discrepancy may actually be due to the fact that your are assuming different distributions for z in the data-generating process and in the model. You use left-truncated log-normal for data generation, whereas you assume that z = m + exp(z_aux), where z_aux follows an improper uniform prior. I'm not sure those two things are the same. Why not model z directly as z ~ lognormal(0,1) [m,]; ?

[Tamas Papp]

Hi,

This is a simplified example from a bivariate censoring problem which
does not admit such a simple formulation, hence the need to model it
this way (see [1]. For the univariate case I could, of course, use
Stan's truncated distros, but the original problem came up for a more
complicated model and I want to keep things similar.

In any case, z_aux has an improper/flat \propto 1 prior, which means
that it does not contribute to the posterior, so it is harmless. I have
a prior on z anyway.

Best,

Tamas

[1] https://groups.google.com/forum/#!topic/stan-users/PZE0tkwJ2A4

On Tue, Aug 19 2014, Luc Coffeng <lucco...@gmail.com> wrote:

> Sorry, I spoke/wrote too soon. Should have read your model code more
> carefully.
>

> The discrepancy may actually may due to the fact that your are assuming

[Luc Coffeng]

I fiddled a little with your model, and realized that the marginal posteriors of m and z_aux are unidentifiable as they added up and both have improper priors. The posterior itself may even be improper because of that (correct me if I'm wrong).

I managed to get the model running using truncated distributions (and supplying initial values, see attachment), but found that after only 500 warm-up simulations, the posterior estimates were sensitive to the actual initial values supplied (R_hat also >> 1). Also, the posterior of m seems to have an implausibly narrow credible interval, which worried me. To quickly check, I cut the number of observations down to 1K (instead of 10K) and cranked up the number of iterations (10K instead of 1K). Interestingly, the model still shows the same behaviour: no convergence, high sensitivity to initial value for m, and still very narrow credible interval for m (maybe this is correct though). Just try setting the inital value for m to 2, and you'll see how the model does not converge within 5000K warm-up samples.

Best,

Luc

[Tamas Papp]

Hi Luc,

Thanks for the work, I really appreciate it! I think that the model does
not have an improper posterior, but just to be on the safe side, I
included a prior on m.

Modifying the code you sent me with a lower bound on z, the model runs
beautifully, quick convergence, etc. However, if I want to do the
transformation manually, I get really bad chains. Examples are attached
(incl. warmup stage).

Either I miscoded something, or Stan is doing something extra that I
didn't include, besided the Jacobian correction. Maybe Bob can help us
out here.

I attach the self-contained code and the plot for m. Any hints on how to
resolve the issue would be appreciated.

(For those not reading the whole thread, I need to transform variables
myself (instead of relying on truncated distros) in a more complex model
with a bivariate censoring region that does not have a rectangular
shape, I ran into problems, hence the univariate exploration.)

Best,

Tamas

[Bob Carpenter]

Stan doesn't do anything extra to the log probability
other than add all the log Jacobian corrections for
constrained parameters. That's it. The rest is up to what's
defined in the model block.

Some well-defined models have the kinds of behaviors you're
seeing in trace plots where parameters get locked into narrow
values. You can try upping the acceptance rate parameter (delta
in all interfaces).

In the by-hand case, I think you're missing the sampling statement:

z[i] ~ lognormal(0,1) T[m, ];


- Bob

> On Tue, Aug 19 2014, Luc Coffeng <lucco...@gmail.com> wrote:
>

> <m-draws.png><stan-censoring-comparison.r>--
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Tamas Papp]

Hi Bob,

It turns out that Stan _does_ do something else, which, of course, needs
to be done: it adjusts for the truncation! The key line that was missing
from my code is

increment_log_prob(-lognormal_ccdf_log(m, 0, 1) * N);

without which the model has the wrong likelihood (cf section 7.6 in BDA,
2nd edition). In the by-hand case, I did have the sampling statement

z ~ lognormal(0,1);

but, of course, compared to

z[i] ~ lognormal(0,1) T[m, ];

the line above was missing. This is an idiotic mistake that I made and
only realized when I wrote out the likelihood explicitly -- this solves
my immediate problem, and I am posting this for the archives. The whole
code is available at

https://gist.github.com/tpapp/fadf244cb8c93244412f

Best,

Tamas

[Bob Carpenter]

Thanks for posting a final solution.

I meant there's nothing else done implicitly to the log probabilities.
The truncations get dropped if the distribution's parameters are constant
and any truncation bounds are constant; otherwise, they affect the
log probability in a way that depends on the parameters.

- Bob