Different results with poisson_log() and poisson_log_glm

Modeling
1

I am using rstan (downloaded this week) to fit a toy data set with 200 rows, three predictors and a response generated from R’s rpois function. I gave the predictors coefficients of 0.3, -0.15, and 0.2 and there is an intercept of 1.4 I fit the data with the model below using either poisson_log or poisson_log_glm. The poisson_log version reports a Total time of 0.5s, Rhats of 1.000 and neff of 2000 to 4000. The version using poisson_log_glm reports a Total time of 11s, Rhats > 140, neff of 2 and a warning that “There were 3997 transitions after warmup that exceeded the maximum treedepth.” The treedepth limit was set to 10. Increasing the max_treedepth to 15 results in a similar warning after a longer run time. The treedepth in the poisson_log version does not exceed 5.

I strongly suspect I am doing something very silly but I cannot see it. Any ideas why the performance is so different?

In case it matters, my R version and system are:
R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux Mint 18.3

data {
    int N; //number obs
    int K; //number of columns in the model matrix
    int y[N]; //response
    matrix[N,K] X; 
  }

parameters {
  real alpha;
  vector[K] beta; //the regression parameters
}
model {  
  alpha ~ normal(0,1.5);
  
  beta ~ normal(0,1);
  
  //y ~ poisson_log(alpha + X * beta); //version 1
  y ~ poisson_log_glm(X, alpha, beta); //version 2
}
2 Likes
2

I’ve had a similar problem. I resolved it by directly updating the target

target += poisson_log_glm_lpmf(y | X, alpha, beta);

rather than using the sampling statement. I think the sampling statement has a bug.

2 Likes
3

@bgoodri is this a RStan-specific issue or is it a problem with math?

4

@aleshing - Yes, that fixed the problem. Thank you! I searched for an issue on github but did not find one. Does one exist?

5

This is the first time I see this reported. Most likely math and not RStan problem. Please make an issue (I’ll be offline next 14 hours)

6

I opened issue #722 in the rstan section.

2 Likes
7

Make one if you don’t see one. If it’s an accidental duplicate, it’s fine, we can clean it up. Post the models and such.

I think the error is here: https://github.com/stan-dev/math/blob/develop/stan/math/prim/mat/prob/poisson_log_glm_lpmf.hpp#L109

include_summand<propto, T_partials_return>::value should be include_summand<propto, T_x, T_alpha, T_beta>::value

The two different specifications compile to cpp differently. The target+ syntax calls:

lp_accum__.add(poisson_log_glm_lpmf<propto__>(y, X, alpha, beta));

and the ~ syntax calls

lp_accum__.add(poisson_log_glm_lpmf<false>(y, X, alpha, beta));

Not sure if this effects any of the other glms.

1 Like
8

Just adding here that this has been confirmed as a bug

9

Thanks to @bbbales2 and @rok_cesnovar, this has been fixed in math and will appear in the upcoming release. :)

2 Likes