Poisson regression model very slow - because of simplex groups of regressors?

[Julian King]

Hi all,

I have a relatively simple Poisson regression problem with n=700, and about 20 regressors (so around 20 regression coefficients).

The adaptation period is, I think, unexpectedly slow, at around one sample per few seconds..

If I start with inits=0,  the chain basically doesn't move for the first few hundred iterations, then moves thereafter.  If I start with init = the parameters from the mode found by the optimiser, the rate of sampling is relatively constant.

After looking at the sampler parameters, max_treedepth is hit the whole way through the first hundreds of iterations.  Increasing max_treedepth up to 12 in accordance with the suggestion in the manual just makes the sampling take a lot longer (tree depth keeps hitting 12).

After adaptation for n=500 iterations (250 warmup and 250 sampling), the effective sample sizes for the parameters are quite high however - often near 250.

I suspect that the problem may be due to some of my regressors.  Groups of regressors unavoidably define a simplex (ie the sum of each of the regressors in a particular group is one).  In fact, I omit one of the regressors from the simplex due to the linear dependence (eg if regressors 1 through 8 define a simplex, then I only use 1 through 7 in the model in order to avoid the degeneracy).  For example, one of my regressors is the percentage of people in a particular age group, eg ages 0-17, 18-29, 30-39, 40-49, 50-59, 60-69, 70-79, and 80+.  In this context I omit the "80+" group to avoid the total degeneracy.

In order to mitigate the potential degeneracy issue and to implement shrinkage, I have put a hierarchical prior on each of the simplex-like groups of regressors with a zero mean (eg age[i] ~ normal(0, sig_age)  in the example above), and I have put a vaguely informative prior on sig_age.

There are a number of possible other causes:

• Scale of regressors (and therefore potential scale of regression coefficients) does vary somewhat - by about 3 orders of magnitude
• I have not used the Matt trick - but if I try removing the hierarchical priors (sig_age in the example above) the speed is not significantly affected, so perhaps this is not the cause

Things I have tried:

• I have vectorised all the sampling statements, so this isn't the issue.
  
• I have tried using a dense mass matrix, and this doesn't seem to make any significant difference.

My question is:  has anyone encountered behaviour like this, and do you have any solutions?  I cannot avoid the use of simplex regressors in this case, but in some cases could amalgamate the groups (but do not wish to do so if I can avoid it).

I appreciate that maybe I just have to live with the slow sampling - perhaps Stan is behaving as expected.  

Thoughts appreciated about where my efforts in optimising the model are best spent.

Thanks in advance!

[Andrew Gelman]

Hi, just a quick guess:  there is typically overdispersion and Poisson regression can run slowly then becuase it’s trying too hard to fit the data.  You could try overdispersed.

A

--
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.

[Bob Carpenter]

Sometimes fitting a more general model can be faster. I'd be surprised if
that were the case with the Poisson --- usually we have more trouble
estimating the overdispersion.

Scale variation of several orders of magnitude essentially introduces
stiffness into Stan's ODE solver --- it can't go too fast in the most
highly constrained dimension.

The problem is that it needs to learn which dimension that is during
warmup, which would explain why things are OK after warmup. If there's a
way you could rescale, that'd probably help a lot.

You can take a look at the scaling parameters after the run to see what
warmup did.

Tree depth of 10 is 1024 log probability and gradient evals per iteration,
and depth 12 is 4 times that many. That's the main reason it's going to
be slow.

- Bob

[Julian King]

Thanks for the suggestions.

I have tried negative binomial with no extra luck.

I will have a look at the sampler parameters after adaptation but the speed up isn't exactly enormous.

I will try one of the continuous approximations to the Poisson and see if that makes a difference. I will try variable rescaling as well.

As an aside, how is the log factorial calculated in the Poisson PDF? Any reason to think there might be a numerical instability there? I guess this wouldn't matter for what I'm doing though as effect should be constant at all iterations

[Bob Carpenter]

To unask the last question, the lgamma function is used for
the Poisson, but it's only calculated if you call poisson_log(...)
as a function. If you use sampling notation, it will be a constant
and be dropped.

When log factorials are needed, we calculate them using the lgamma function
from Boost:

http://www.boost.org/doc/libs/1_59_0/libs/math/doc/html/math_toolkit/sf_gamma/lgamma.html

and digamma and trigamma for derivatives.

- Bob

[Julian King]

Thanks Bob - much appreciated.

The suggestion to look at the mass matrix was a good one. The inverse diagonal weights reported range from ~1e-6 to ~10. Looking at the posterior summary outputs, it turns out that some of my parameters are just much better constrained in the regression than others.

I guess that this is the source of the problem. I will try rescaling my inputs so that the likely posterior has a much more uniform scale over the different variables.

Thanks once again for your help guys - I realise how busy you are and how many requests you get and am grateful for your assistance.

Is there any way to pass in an initial guess at the mass matrix? Or is this for a future release? Being able to specify guessed scales for the parameters would be more convenient in some applications than having to recode parameters.

[Julian King]

[Julian King]

[Bob Carpenter]

> On Sep 26, 2015, at 8:50 PM, Julian King <jking...@gmail.com> wrote:
>
> Thanks Bob - much appreciated.
>
> The suggestion to look at the mass matrix was a good one. The inverse diagonal weights reported range from ~1e-6 to ~10. Looking at the posterior summary outputs, it turns out that some of my parameters are just much better constrained in the regression than others.
>
> I guess that this is the source of the problem. I will try rescaling my inputs so that the likely posterior has a much more uniform scale over the different variables.

That should help considerably.

> Thanks once again for your help guys - I realise how busy you are and how many requests you get and am grateful for your assistance.

You're welcome. Problems like this one are good because
we need to diagnose how Stan can fail. Ideally, we'd like to
improve our adaptation routines as well as our diagnostics.

> Is there any way to pass in an initial guess at the mass matrix? Or is this for a future release? Being able to specify guessed scales for the parameters would be more convenient in some applications than having to recode parameters.

Not right now, but there will be in the future. It's
part of the Stan 3 refactor of the API between the
interfaces (RStan, PyStan, CmdStan) and the language.
I'm not forgetting the other interfaces, but they all
call CmdStan, so should be insulated from these changes.

- Bob

[Julian King]

An update on this:

• I tried a negative binomial model, then just moved back to plain Poisson because the overdispersion in the data is not material
• Placing somewhat more informative (but still pretty vague) priors seems to have helped
• I am using a dense mass matrix and this helps, I think
• I have used adapt_delta = 0.9, and this seems to help
  
• Looking at the posterior for coefficients of the simplex regressor groups reveals very high correlations (r > 0.9), so this isn't exactly a very easy posterior to sample from.  Somewhat informative priors seem to have helped here, as per the suggestion in the manual for nearly-unidentified parameters.
  
• It seems that simply a lot of warmup iterations are needed in order for Stan to learn the shape of the posterior.  With n=500 warmup iterations, max_treedepth is hit on most iterations early on, and progressively less later.  With larger numbers of warmup iterations the average number of leapfrog steps seems to decline, but there doesn't appear to be a 'eureka' moment where suddenly sampling quality drastically improves.

Unfortunately it just looks like I have a slow model.  My data re are constructed so that y=0 never appears, so I have to use a ~ poisson T[1,] statement, which is unfortunately not vectorised.

I'm happy that the model is working as intended, and the run time is long but not unmanageably.   Thanks for your help.

[Bob Carpenter]

I don't know how strict that truncation is (seems to be Poisson and
overdispersed and zero-inflated Poisson week), but rather than

y[n] ~ poisson(alpha[n]) T[1,];

you could use

(y[n] - 1) ~ poisson(alpha[n]);

Of course, the interpretation of alpha[n] is slightly different.

Alternatively, you can do the simple case of the Poisson by realizing
that the relevant CDF (1 - poisson_cdf(1, alpha[n])) can be coded
as (1 - exp(poisson_log(0, alpha[n])). Because you're going to need
that on the log scale, and it's a common operation, we have it built in
as

log1m_exp(poisson_log(0, alpha[n]);

So you could code the whole model as

y ~ poisson(alpha);
for (n in 1:N)
increment_log_prob(log1m_exp(poisson_log(0, alpha[n]));

Considerably uglier. Now if the alpha[n] were all the same, it
could be made much more efficient. As is, allowing the vectorization
now in the Poisson sampling statement shouldn't be too bad.

- Bob

[Julian King]

Shouldn't this be

y ~ poisson(alpha);
for(i in 1:n){
    increment_log_prob(-log1m_exp(-alpha[i]));
}
??

According to p109 of the manual, the PDF should be
p(y | lambda) = Poisson(y | lambda) / (1 - PoissonCDF(0 | lambda))
so log p = log_Poisson(y | lambda) - log(1 - PoissonCDF(0 | lambda))
 = log_Poisson(y | lambda) - log(1 - exp(-lambda)
 = log_Poisson(y | lambda) - log1m_exp(-lambda)

[Andrew Gelman]

Hi, once this is all working it would be great to write up as a short vignette, as this general problem arises a lot. We can publish it in the Stan Journal.
A

[Bob Carpenter]

??? There's nothing in the brackets. Did something get lost in translation?
You don't need brackets around a single statement if that's what you're worried
about.

I think what I wrote below is OK assuming alpha's computed ahead of time
as a vector.

- Bob

> On Sep 29, 2015, at 7:42 AM, Julian King <jking...@gmail.com> wrote:
>
> Shouldn't this be
>
> y ~ poisson(alpha);
> for(i in 1:n){
>
> }
>

[Julian King]

Bob,

Sorry, I accidentally hit enter on Google Groups and then edited it after.  I assume you're getting this through the mailing list, so didn't see the edit.  What I wrote was this:
-----

y ~ poisson(alpha);
for(i in 1:n){

    increment_log_prob(-log1m_exp(-alpha[i]));

}
??

According to p109 of the manual, the PDF should be
p(y | lambda) = Poisson(y | lambda) / (1 - PoissonCDF(0 | lambda))
so log p = log_Poisson(y | lambda) - log(1 - PoissonCDF(0 | lambda))
 = log_Poisson(y | lambda) - log(1 - exp(-lambda)
 = log_Poisson(y | lambda) - log1m_exp(-lambda)

----

[Bob Carpenter]

Clever --- I hadn't thought to reduce poisson_log(0, lambda)
further. It inverts the "zeroes trick" from JAGS. :-)

I'll add a note that you can keep going on the reduction to the manual
and thank you for pointing it out:

https://github.com/stan-dev/stan/issues/1617#issuecomment-144304649

- Bob

[Julian King]

No problems.  Not sure how it affects the vectorisation and so whether there is any significant speed difference.  Obviously if lambda is the same everywhere you can make it common, but in a Poisson regression model it's not.

[Julian King]

Just to provide a further update on my general problem:

• I have not been able to materially speed up the model (iterations per second) through simply recoding the Poisson calculation
• I continued to struggle with convergence issues, until I standardised every regressor to have mean zero and standard deviation 1.  This has made a massive improvement to the model, I suspect by reducing posterior correlations (I have interaction terms as well, which behave oddly unless standardised).  I had some concerns about this given the simplex-like nature of blocks of my regressors, but in fact this seems to be fine.
• I have had to use adapt_delta = 0.9 in order to minimise the number of divergent transitions
• I have increased max_treedepth to 11, but still hit up against the max_treedepth quite often. 
  

[Michael Betancourt]

Have you plotting the parameter correlations? This sounds a lot like non-identifiability (your posterior is concentrating

on a warped cigar, which requires very small step sizes and long integration times to explore) which isn’t surprising given 

the simplex constraint.

[Julian King]

Yes, the parameter correlations are very high given the simplex constraint.  I have a Gaussian shrinkage prior on every parameter to try to deal with the problem, so the prior is nominally proper, but appreciate that it can take a long time to explore.

With the shrinkage prior, the posterior is not entirely elliptical, which adds to the problem.

Warmup of 1,000 to 2,000 iterations seems to produce reasonable results thereafter (although with about one effective sample per 10-20 iterations).  It's not slow, but thankfully I can run a large number of chains in parallel, so it's not intolerable to get a good level of accuracy.  It makes debugging a pain though.

Thanks for the comment.

[Bob Carpenter]

Might it be possible to fix this by going to a (K-1)-factor
basis for the simplex rather than a K-factor basis? We need more
discussion of this issue in the manual. Our K-simplex parameters
in Stan are parameterized in the unconstrained space with (K-1) free
parameters.

- Bob

[Julian King]

Bob,

For each of my groups of regressors which sum to one, I omit one of them in order to deal with this issue.  Eg if I have 7 regressors in a group which sum to one, I only use 6 of them.

Note that the parameters are not constrained to sum to one, only the regressors.

Julian

[Julian King]

Another update:  I have tried non-centred parameterisations for the simplex-regressor parameters.  This vastly improves the mixing on these, but shifts the problem to the regularisation priors that I impose on them, as well as other parameters in the model (these mix poorly initially, and apparently require long run times in order to learn the relevant scales).