Iteration confusion with zero inflated poisson model

[Timothy Mastny]

Hello,

I am trying to fit a zero inflated poisson model, very similar to the one in the distributional modes vignette: https://cran.r-project.org/web/packages/brms/vignettes/brms_distreg.html

The main difference is that I am trying to include an exposure/offset. Fitting the model without the offset goes quickly with no trouble (low Rhat).

At 2000 iterations, the model appears to fit fine, with Rhat=1 and sensible parameter estimates. When I try 4000 iterations, the model has high Rhat values and says that max_treedepth has been reached. I tried the model on the default tree depth of 10 and 12. I tried with 15 but after 45 minutes the model wouldn't fit.

Furthermore, fitting the model in the rethinking package, even with 4000 iterations, I get similar parameter estimates as the 2000 iteration model. 

Anyone have any perspective what might be going on? Have a misspecified the brms model somehow? Any advice on debugging?

Here is the minimum working code and output to show this effect:

library(brms)
library(rethinking) # for data
data(Fish)
d <- Fish


rstan_options (auto_write=TRUE)
options (mc.cores=parallel::detectCores ()) # Run on multiple cores


mod.brms.2000 <- brm(bf(fish_caught ~ persons + camper + child + offset(log(hours)),
                   zi ~ child),
                data=d, family=zero_inflated_poisson(),
                prior = c(set_prior("normal(0,10)", class = 'b'),
                          set_prior("normal(0,10)", class = 'Intercept')),
                iter=2000)


summary(mod.brms.2000)


#  Family: zero_inflated_poisson
#   Links: mu = log; zi = logit
# Formula: fish_caught ~ persons + camper + child + offset(log(hours))
#          zi ~ child
#    Data: d (Number of observations: 250)
# Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#          total post-warmup samples = 4000
#     ICs: LOO = NA; WAIC = NA; R2 = NA
#
# Population-Level Effects:
#              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
# Intercept       -2.20      0.17    -2.53    -1.89       2967 1.00
# zi_Intercept    -1.42      0.32    -2.10    -0.85       4000 1.00
# persons          0.74      0.04     0.65     0.83       2852 1.00
# camper          -0.39      0.09    -0.57    -0.20       3218 1.00
# child            0.59      0.09     0.40     0.77       2245 1.00
# zi_child         0.71      0.46    -0.39     1.52       1898 1.00
#
# Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
# is a crude measure of effective sample size, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).


mod.brms.4000 <- brm(bf(fish_caught ~ persons + camper + child + offset(log(hours)),
                   zi ~ child),
                data=d, family=zero_inflated_poisson(),
                prior = c(set_prior("normal(0,10)", class = 'b'),
                          set_prior("normal(0,10)", class = 'Intercept')),
                iter=4000, control=list(max_treedepth=12))


summary(mod.brms.4000)


# Warning messages:
#   1: There were 2385 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 12. See
# http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
# 2: Examine the pairs() plot to diagnose sampling problems




#  Family: zero_inflated_poisson
#   Links: mu = log; zi = logit
# Formula: fish_caught ~ persons + camper + child + offset(log(hours))
#          zi ~ child
#    Data: d (Number of observations: 250)
# Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#          total post-warmup samples = 8000
#     ICs: LOO = NA; WAIC = NA; R2 = NA
#
# Population-Level Effects:
#              Estimate Est.Error  l-95% CI u-95% CI Eff.Sample Rhat
# Intercept       -2.24      0.17     -2.59    -1.93        115 1.06
# zi_Intercept    -1.36      0.31     -2.01    -0.81        281 1.02
# persons          0.75      0.05      0.66     0.84         96 1.06
# camper          -0.38      0.10     -0.56    -0.19        148 1.01
# child            0.48      0.14      0.23     0.74          6 1.54
# zi_child     -2985.21   4473.48 -14226.41     1.35          7 1.97
#
# Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
# is a crude measure of effective sample size, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).
# Warning message:
#   The model has not converged (some Rhats are > 1.1). Do not analyse the results!
#   We recommend running more iterations and/or setting stronger priors.


d$loghours <- log(d$hours)
m11h6.baseline <- map2stan(
  alist(
    fish_caught ~ dzipois(p, lambda),
    logit(p) <- zi + zC*child,
    log(lambda) <- loghours + al + bC*camper + bP*persons + bCh*child,
    c(zi, al) ~ dnorm(0, 10),
    c(zC, bC, bP, bCh) ~ dnorm(0, 10)
  ), data = d, chains = 4, warmup = 1000, iter = 4000)


precis(m11h6.baseline)


##      Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## zi  -1.42   0.31      -1.89      -0.92  7835 1.00
## al  -2.20   0.16      -2.44      -1.93  6009 1.00
## zC   0.61   1.01      -0.04       1.42   551 1.01
## bC  -0.39   0.09      -0.54      -0.24  8645 1.00
## bP   0.74   0.04       0.67       0.81  6261 1.00
## bCh  0.59   0.10       0.44       0.74  4913 1.00

Thanks,

Tim

[Paul Buerkner]

Looks like this may be bad luck with your chain. Try changing argument seed of brm.

--
You received this message because you are subscribed to the Google Groups "brms-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+unsubscribe@googlegroups.com.
To post to this group, send email to brms-...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/be0c3229-8e4d-476c-87f7-0b34f125e14e%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Timothy Mastny]

I tried fitting the model a few different times using the random seed method described in the brms manual. The mean value of zi_child is less extreme, but still has a very large Rhat.

I'm also still getting this warning about max tree depth:

Warning messages:

1: There were 2385 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 12. See

http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

2: Examine the pairs() plot to diagnose sampling problems

Stan also recommends examining the pairs plot. The marginal posterior distribution for b_zi_child does look very weird. The weirdness disappears when I remove the offset term, which makes me think there is something going on with that. Here is the output without the offset:

mod.brms <- brm(bf(fish_caught ~ persons + camper + child, 
                   

zi ~ child), 
                data=d, family=zero_inflated_poisson(),
                prior = c(set_prior("normal(0,10)", class = 'b'),
                          set_prior("normal(0,10)", class = 'Intercept'))

,
                iter=4000, seed = sample(1e+7, size = 1))

summary(mod.brms)

##  Family: zero_inflated_poisson 
##   Links: mu = log; zi = logit 
## Formula: fish_caught ~ persons + camper + child 
##          zi ~ child
##    Data: d (Number of observations: 250) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1; 
##          total post-warmup samples = 8000
##     ICs: LOO = NA; WAIC = NA; R2 = NA
##  
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept       -1.07      0.18    -1.44    -0.72       6179 1.00
## zi_Intercept    -0.95      0.26    -1.52    -0.48       7122 1.00
## persons          0.89      0.05     0.80     0.98       6023 1.00
## camper           0.77      0.09     0.59     0.96       7394 1.00
## child           -1.17      0.09    -1.36    -0.99       6928 1.00
## zi_child         1.22      0.27     0.70     1.77       8000 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Here is an example of one of the runs with a new random seed:

mod.brms <- brm(bf(fish_caught ~ persons + camper + child + offset(log(hours)), 
                   

zi ~ child), 
                data=d, family=zero_inflated_poisson(),
                prior = c(set_prior("normal(0,10)", class = 'b'),
                          set_prior("normal(0,10)", class = 'Intercept'))

,
                iter=4000, seed = sample(1e+7, size = 1))

summary(mod.brms)

##  Family: zero_inflated_poisson 
##   Links: mu = log; zi = logit 
## Formula: fish_caught ~ persons + camper + child + offset(log(hours)) 
##          zi ~ child
##    Data: d (Number of observations: 250) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1; 
##          total post-warmup samples = 8000
##     ICs: LOO = NA; WAIC = NA; R2 = NA
##  
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept       -2.23      0.17    -2.56    -1.92         70 1.05
## zi_Intercept    -1.39      0.30    -2.03    -0.84       8000 1.01
## persons          0.75      0.04     0.66     0.83         15 1.07
## camper          -0.39      0.09    -0.57    -0.21        306 1.01
## child            0.54      0.13     0.24     0.76          4 1.44
## zi_child      -533.02   1080.18 -3609.19     1.42          3 2.51
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+...@googlegroups.com.

[Paul Buerkner]

Try setting proper priors in the zi coefficients for instance via

prior = prior(normal(0, 10), dpar = zi) + prior(normal(0, 10), class = Intercept, dpar = zi)

To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+unsubscribe@googlegroups.com.

To post to this group, send email to brms-...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/8b57cc7d-3087-45bd-9cef-a263803f43e8%40googlegroups.com.

[Timothy Mastny]

Thanks that was my issue. It fit right away.

For future readers, my problem was that I thought I had put the normal priors on the zi coefficients. I wasn't aware of the distribution parameter argument dpar.

Thanks for the help,

Tim