Zero-truncated lognormal-Poisson

[Janne Sinkkonen]

I have trouble with this simple model:

data {

  int N;

  int freqs[N];

}

parameters {

  real mu;

  real<lower=0> sigma;

  vector<lower=0>[N] theta;

  

}

model {

  theta ~ lognormal(mu, sigma);

  freqs ~ poisson(theta);

  for (n in 1:N) increment_log_prob(-log1m_exp(-theta[n]));

  //  mu ~ normal(-1, 6);

  //  sigma ~ lognormal(1, 6);

}

It is supposed to be a compound lognormal-Poisson model, for frequency data (freqs) with zero frequencies not observed.

So data typically looks like this:

    80 69 39 31 29  9  7  7  6  5  4  4  3  3  3  2  2  2  1  1  1  1  1  1  1  1  1

.. but is longer, so that its length (N) is around 1000-2000. 

I see a systematic bias with real data if I make a posterior check by generating draws from (mu, sigma). Also the estimates are clearly off for simulated data, and different from ML estimates from R package sads.

The problem seems to be the censoring and the likelihood correction, as I tried simulated data with vanilla lognormal-Poisson and got good estimates. I have tried to write the likelihood correction with poisson_ccdf as well, and directly without log1m_exp(), with no change.

There is probably something obviously wrong with the model or the approach to correct for missing zeros that I don't see, but what is it? :)

(An equivalent censored gamma-Poisson model with theta integrated out seemed to work fine, although I didn't test it very thoroughly. My Stan is from 2.9.0 from CRAN, installed today, and the compiler is clang.)

[Bob Carpenter]

That truncation calculation looks right to me.

log(1 / (1 - poisson_cdf(0 | theta[n])))
= -log(1 - poisson_cdf(0 | theta[n]))
= -log(1 - poisson(0 | theta[n]))
= -log(1 - exp(-theta[n]))
= -log1m_exp(-theta[n])

What kind of systematic bias are you seeing? Are
the intervals right? Does the limit as you get
more data do the right thing? Are multiple chains converging
OK as measured by r-hat and are you getting reasonable n_eff
estimates?

- Bob

> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Janne Sinkkonen]

Comments below.

On Thursday, June 16, 2016 at 7:39:37 PM UTC+3, Bob Carpenter wrote:

What kind of systematic bias are you seeing?

Like with data generated from mu=-1, sigma=3, I get estimates mu≈1.37±0.05 (sd)  and sigma=2.06±0.04, while ML from the sads package are for example -0.60±0.25 and 2.78±0.11 (My impression is that even the ML estimates are off the ±2sd limits more than 5% of the time, but they are on the right ballpark at least, and better with larger data.)

Compared to the ML estimates, Intervals from the Stan model seem overly tight. 

Making data smaller in the sense of number of counts (N) leads to wider intervals and maybe a bit better estimates, but the intervals are still maybe too narrow. (And mixing gets poor, so it's hard to see what's happening near the prior.)

Increasing data tightens the intervals, but keeps the estimates off (pretty much where they were, so with increasing N it converges into wrong estimates with overly tight intervals). 

Making the true sigma smaller makes the estimates closer to true values (which is understandable for the lognormal then kind of vanishes), but mixing gets worse so I can't take the sigma very small. I haven't played with the true mean much.

Posterior means of theta are probably ok, with a downward bias when freq is one or small. (This is ok I guess, for then the lognormal prior gives more evidence than data. At larger values of freq, data totally overrides the prior as it should with Poisson.)

(It somehow it feels like the lognormal is the culprit, but I'm mostly just baffled...)

Chains are ok as far as I can see, in terms of trace plots, r-hat and n_eff. A few (≈ 1–5) numerical problems or divergent transitions are reported. 

I tested with bare Poisson, the theta removed but censoring left in, and then the estimates are right. I even tested a lognormal with data and the Poisson thing left out, and it also seems to be fine. :) 

All this is with my Mac, but there's an AWS instance available where I could try to run the model later. 

Files used for tests are attached below. 

[Bob Carpenter]

You could try translating poisson-lognormal to a regular
normal scale and transforming the result. That'd at least
verify that all the coding is right.

In your model, it'd be

parameters {
real mu;
real log_sigma;
real log_theta[N];
}
model {
log_theta ~ normal(mu, exp(log_sigma));
freqs ~ poisson_log(log_theta);

for (n in 1:N) increment_log_prob(-log1m_exp(-theta[n]));

mu ~ normal(-1, 6);
log_sigma ~ normal(3, 6);
}


That should be the same model, only ith the parameters now expressed
on the unconstrained scale.

And to double check the increment_log_prob version, I'd try

for (n in 1:N)

freqs[n] ~ poisson_log(log_theta) T[1, ];

I used poisson_log(log_theta), which is equivalent to
poisson(exp(log_theta)).

Are those posterior sd or MCMC std error? We usually try
to calibrate with 80% or 90% intervals, or even 50% intervals
if there are enough values.

I always get confused on where the logs are happening in the
lognormal, so I'd check that; the inner term for the
lognormal log density is

((log y - mu) / sigma)^2

- Bob

> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <poisson-lognormal-debug.R><poisson-lognormal.stan><poisson.stan>

[Janne Sinkkonen]

Change to normal() and log_sigma did not change anything. The likelihood correction term needed an extra exp(),

 for (n in 1:N) increment_log_prob(-log1m_exp(exp(-log_theta[n]))); 

that was not in your changes. 

The truncation syntax then: it destroys any mixing of the chains (5000 iterations):

(I have actually tried truncation syntax before, but didn't mention it, because it is a bit separate issue.) Does this mean that the likelihood correction is wrong? Well, the same likelihood correction term works with a bare Poisson model without the normal/lognormal layer. I haven't tried truncation syntax there. 

BTW, these are the chains with the manual correction term (500 iterations):

My sd's are posterior standard deviations, and they were there to compare with the ML estimate which gives approximate sd's. The (false) posteriors seem normal enough to make the sd estimates robust. 

[Janne Sinkkonen]

I don't see the figs in my previous mail, so post them here again as attachments. 

[Janne Sinkkonen]

Not my day with attachments... here is the other one with manual truncation term and good mixing.