poisson_log_glm sampling statement fails to fit

[FJCC]

Summary:

Fitting data with
y ~ poisson_log_glm(X, alpha, beta)
fails with almost 100% treedepth warnings but using
target += poisson_log_glm_lpmf(y | X, alpha, beta)
works well

Description:

Poisson regression on a toy data set using
y ~ poisson_log_glm(X, alpha, beta)
in the model section returns a warning for 3997 (out of 4000) transitions exceeding the maximum treedepth, Rhat is > 100 and neff = 2. Using either
y ~ poisson_log(alpha + X * beta)
or
target += poisson_log_glm_lpmf(y | X, alpha, beta)
results in a high quality fit: Rhat ~ 1.000, neff in the range 2000 -4000. The rstan code for generating the toy data and fitting with the three versions is below.

Reproducible Steps:

set.seed(1)
DF2 <- data.frame(X1 = rnorm(200, 0, 1), X2 = runif(200, -2, 2), X3 = rnorm(200, 1, 2))
DF2 <- DF2 %>% mutate(Y = 0.3 * X1 - 0.15 * X2 + 0.2 * X3 + 1.4)
DF2 <- DF2 %>% mutate(Cnt = rpois(200, lambda =  exp(Y)))

dataList2 <- list(N = nrow(DF2), K = 3, y = DF2$Cnt, X = as.matrix(DF2[, 1:3]))
PoissonMod <- "
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, works well
  //y ~ poisson_log_glm(X, alpha, beta); //version 2, treedepth failures near 100%
  target += poisson_log_glm_lpmf(y | X, alpha, beta); //version 3, works well
}
"
PoisOut <- stan(model_code = PoissonMod, data = dataList2, verbose = FALSE, control = list(max_treedepth = 10))

RStan Version:

rstan_2.19.2 StanHeaders_2.19.0

R Version:

R version 3.5.2 (2018-12-20)

Operating System:

Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux Mint 18.3

[bob-carpenter]

Thanks for reporting and providing a careful comparison. That indeed looks like a bug somewhere.

[bbbales2]

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));

I don't know if the difference in the calls is a bug or not. Someone else can take this. Just this is what I found looking earlier.

[avehtari]

Pinging @t4c1 and @rok-cesnovar who have recently been working on compound functions

[rok-cesnovar]

For now, I can confirm that I can reproduce this on the current develop of cmdstan which means that its still an issue in the latest Stan Math. I went to double-check because there was a lot of churn on GLMs since Stan Math 2.19 which you are using via Stan Headers 2.19.

I will explore this further but the code to reproduce with cmdstanr is this (put the stan code in a separate file:

library(dplyr)
library(cmdstanr)

setwd("~/Desktop/Stan/tests/")
set_cmdstan_path("/home/rok/Desktop/Stan/cmdstan/")
file <- file.path("poisson_glm_bug.stan")
mod <- cmdstan_model(file)

set.seed(1)
DF2 <- data.frame(X1 = rnorm(200, 0, 1), X2 = runif(200, -2, 2), X3 = rnorm(200, 1, 2))
DF2 <- DF2 %>% mutate(Y = 0.3 * X1 - 0.15 * X2 + 0.2 * X3 + 1.4)
DF2 <- DF2 %>% mutate(Cnt = rpois(200, lambda =  exp(Y)))
dataList2 <- list(N = nrow(DF2), K = 3, y = DF2$Cnt, X = as.matrix(DF2[, 1:3]))

fit <- mod$sample(
   data = dataList2,
   num_chains = 4,
   num_cores = 1,
   max_depth = 10
)

Wit the ~ version I get the same issue:

Total execution time: 15.5 seconds.
3997 of 1000 (400%) transitions hit the maximum treedepth limit of 10 or 2^10-1 leapfrog steps.
Trajectories that are prematurely terminated due to this limit will result in slow exploration.
Increasing the max_depth limit can avoid this at the expense of more computation.
If increasing max_depth does not remove warnings, try to reparameterize the model.

[rok-cesnovar]

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

Can confirm that this fixes it.

[bbbales2]

Actually I think that include_summand doesn't need to be there at all cause this is already before it: https://github.com/stan-dev/math/blob/develop/stan/math/prim/mat/prob/poisson_log_glm_lpmf.hpp#L76

I think the propto thing is another bug as well.

Like, the _lpmf should have 'propto = false' cause in those cases we want the full lpmf, not just something up to a proportionality constant suitable for MCMC (which is what propto = true means). This is important for like the loo calculations and such in generated quantities.

For the beta parameter (which I assume is correct) I get:

lp_accum__.add(normal_lpdf<false>(beta, 0, 1));

for

target += normal_lpdf(beta | 0, 1);

and

lp_accum__.add(normal_log<propto__>(beta, 0, 1));

for

beta ~ normal(0, 1);

[rok-cesnovar]

I have the fix up on a branch together with the same fix for the GPU implementation, just need to write a test that would fail now and pass with the fix. Should be ready tomorrow.

[bob-carpenter]

Like, the _lpmf should have 'propto = false' cause in those cases we want the full lpmf, not just something up to a proportionality constant suitable for MCMC (which is what propto = true means).

That's right---using foo_lpfm(...) in Stan should have propto=false generated for the call.

[rok-cesnovar]

@bbbales2 With cmdstan+stanc3 I get the following:

target += poisson_log_glm_lpmf(y | X, alpha, beta);
produces
lp_accum__.add(poisson_log_glm_lpmf<false>(y, X, alpha, beta));

and

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

On cmdstan 2.21 (with the c++ stan compiler) I get the following:

target += poisson_log_glm_lpmf(y | X, alpha, beta);
produces
lp_accum__.add(poisson_log_glm_log(y, X, alpha, beta));

and

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

So I think this is all fine? But do please confirm. I also opened the PR with the fix for the Math code.