Lack of trace convergence and low n_eff for poisson models with interaction terms

[George Peck]

Hello,

I am fitting count data to various poisson models using the code and packages that integrate with rstan are outlined in chapter 10 of "Statistical Rethinking" by Richard McElreath,  I have constructed models for my count data, and run them in RStudio using rstan with the ancillary functions from McElreath's github site (see http://xcelab.net/rm/software/).   For the first analysis I constructed models using three main factors.  The full model has eight terms: intercept, three main terms, three two-way interactions, and one three way interaction.  The Widely Applicable Information Criterion (WAIC) was highest for the full model (all interaction terms), and second highest for the model with no three-way interaction.  Interestingly, the next best model was the one with only the intercept and the trap term.  However, when I checked the MCMC diagnostics (Markov chain trace plot convergence and the effective number of samples (n_eff; pages 253-255 in McElreath)), there was no trace convergence and very low n_eff for the full model and the model with all terms except the three-way interaction.  McElreath suggests manipulating priors to be "weakly informative", but I am not sure how to do that.

Here is my model code:

#model with all two-way and three-way interactions

m10.10 <- map(

  alist(

    count ~ dpois( lambda ),

    log(lambda) <- a + bt*trap + 

      bd*date + bs*site + btd*trap*date + bts*trap*site +

      bds*date*site + btds*trap*date*site,

    a ~ dnorm(0,100),

    bt ~ dnorm(0,1),

    bd ~ dnorm(0,1),

    bs ~ dnorm(0,1),

    btd ~ dnorm(0,1),

    bts ~ dnorm(0,1),

    bds ~ dnorm(0,1),

    btds ~ dnorm(0,1)

  ),

  data = d)

# Model with no three way interaction

m10.11 <- map(

  alist(

    count ~ dpois( lambda ),

    log(lambda) <- a + bt*trap + 

      bd*date + bs*site + btd*trap*date + bts*trap*site +

      bds*date*site,

    a ~ dnorm(0,100),

    bt ~ dnorm(0,1),

    bd ~ dnorm(0,1),

    bs ~ dnorm(0,1),

    btd ~ dnorm(0,1),

    bts ~ dnorm(0,1),

    bds ~ dnorm(0,1)

  ),

  data = d)

# Model with no date*site interaction

m10.12 <- map(

  alist(

    count ~ dpois( lambda ),

    log(lambda) <- a + bt*trap + 

      bd*date + bs*site + btd*trap*date + bts*trap*site,

    a ~ dnorm(0,100),

    bt ~ dnorm(0,1),

    bd ~ dnorm(0,1),

    bs ~ dnorm(0,1),

    btd ~ dnorm(0,1),

    bts ~ dnorm(0,1)

  ),

  data = d)

# Model with no bts, bds, and btds interactions

m10.13 <- map(

  alist(

    count ~ dpois( lambda ),

    log(lambda) <- a + bt*trap + 

      bd*date + bs*site + btd*trap*date,

    a ~ dnorm(0,100),

    bt ~ dnorm(0,1),

    bd ~ dnorm(0,1),

    bs ~ dnorm(0,1),

    btd ~ dnorm(0,1)

  ),

  data = d)

# Model with no interactions

m10.14 <- map(

  alist(

    count ~ dpois( lambda ),

    log(lambda) <- a + bt*trap + 

      bd*date + bs*site,

    a ~ dnorm(0,100),

    bt ~ dnorm(0,1),

    bd ~ dnorm(0,1),

    bs ~ dnorm(0,1)

  ),

  data = d)

#

#submodel 15 (omit location)

m10.15 <- map(

  alist(

    count ~ dpois( lambda ),

    log(lambda) <- a + bt*trap + bd*date,

    a ~ dnorm( 0 , 100 ),

    c(bt,bd) ~ dnorm( 0 , 1 )

  ), data=d)

#submodel 16 (omit date)

m10.16 <- map(

  alist(

    count ~ dpois( lambda ),

    log(lambda) <- a + bt*trap,

    a ~ dnorm( 0 , 100 ),

    bt ~ dnorm( 0 , 1)

  ), data=d )

#submodel 17 (intercept only)

m10.17 <- map(

  alist(

    count ~ dpois( lambda ),

    log(lambda) <- a,

    a ~ dnorm(0,100)

  ), data=d)

### Here is my check for trace convergence and n_eff:

### Verification that the MAP estimates accurately describe the shape of the posterior distribution

## Checking n effective and trace plots for convergence

m10.10stan <- map2stan( m10.10 , iter=3000 , warmup=1000 , chains=4 )

precis(m10.10stan)

pairs(m10.10stan)

plot(m10.10stan)

# neff very small and traces do not converge

m10.11stan <- map2stan( m10.11 , iter=3000 , warmup=1000 , chains=4 )

precis(m10.11stan)

pairs(m10.11stan)

plot(m10.11stan)

#neff very small and traces do not converge

m10.12stan <- map2stan( m10.12 , iter=3000 , warmup=1000 , chains=4 )

precis(m10.12stan)

pairs(m10.12stan)

plot(m10.12stan)

#neff large, traces converge during warmup

m10.13stan <- map2stan( m10.13, iter=3000 , warmup=1000 , chains=4)

precis(m10.13stan)

pairs(m10.13stan)

plot(m10.13stan)

# neff large, traces converge during warmup

m10.14stan <- map2stan( m10.14, iter=3000 , warmup=1000 , chains=4)

precis(m10.14stan)

pairs(m10.14stan)

plot(m10.14stan)

# neff large, traces converge during warmup

m10.15stan <- map2stan( m10.15, iter=3000 , warmup=1000 , chains=4)

precis(m10.15stan)

pairs(m10.15stan)

plot(m10.15stan)

# neff large, traces converge during warmup

m10.16stan <- map2stan( m10.16 , iter=3000 , warmup=1000 , chains=4 )

precis(m10.16stan)

pairs(m10.16stan)

plot(m10.16stan)

# neff large, traces converge during warmup

m10.17stan <- map2stan( m10.17 , iter=3000 , warmup=1000 , chains=4 )

precis(m10.17stan)

pairs(m10.17stan)

plot(m10.17stan)

# neff large, trace converges during warmup

### Here is the code for WAIC computation and graphing

#Model comparisons

#Compare all submodels using "Widely Applicable Information Criterion" or WAIC

#Adding n=le4 for more stable comparison

#Will also plot the comparison

(traps.compare <- compare(m10.10,m10.11,m10.12,m10.13,m10.14,m10.15,m10.16,m10.17 ,n=1e4))

plot(traps.compare)

###  Thank you for the help!!!

George

[Richard McElreath]

The priors don't seem too bad, although you might try making dnorm(0,100) into dnorm(0,10).

My hunch is the scale on your predictors. If date, for example, has a large range, then the model will have trouble adapting, given that your priors on all the coefficients imply the same scale for all the predictors. Try standardizing your predictors to the same scale.

I also wonder whether "trap" is meant to be an index variable rather than a continuous predictor. But that's a secondary issue.

[George Peck]

Hi Richard,

Making dnorm(0,100) into dnorm(0,10) allowed trace convergence during warmup and n_eff was quite large (thank you).  However, the pairs plot revealed obvious correlation among pairs of parameters. I got this error message:

Warning messages:

1: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. 

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

 

Not sure what "adapt_delta" is.  Will standardizing predictors (centering on zero; R code 10.50 in your text) help remove this error?

While trap is categorical, I am comfortable treating it as continuous, in the sense that trap attractiveness is continuous (my main interest).  Date and location are coded 1,2,3,4 in my matrix, but I could center. 

Please advise.  

[Bob Carpenter]

Whether the divergences constitute an "error" or a "warning" depends
on who you talk to :-) adapt_delta is an underlying Stan control
variable for the target acceptance rate (related to Metropolis
acceptance, but quite different)---the main point is that an adapt_delta
closer to one means the algorithm is more careful in simulating the
Hamiltonian (if the simulation is perfect, the acceptance rate is 1),
which in practice reduces to smaller step sizes because we use the
gradient to follow the posterior in small steps (that's the numerical
integrator to solve the Hamiltonian diff eqs presented by HMC).

We've had a hard time reporting warnings in Stan that users don't
interpret as errors. Is that because most R programs don't report
any warnings at all?

Yes, standardizing predictors can help, but the usual culprit is
the posterior geometry, as in Neal's funnel example (see the Stan
manual for a description).

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

[Richard McElreath]

Great, sounds like all is well now.

It's not necessarily a problem that parameters are correlated in the posterior. That's normal. It's just that chains with high correlations tend to mix less efficiently.

That one divergent transition is from the final step taken to compute DIC. Scroll up and you'll perhaps see that.

Both divergent transitions and the adapt_delta control parameter are addressed in later chapters, I think in Chapter 13 and Chapter 14, if not already in Chapter 12. It's a Stan setting, so if you search this list, you'll find lots of examples of its use.

And yes, I'd standardize (not just center) all continuous predictors, just so the priors are easy to interpret and on the same conceptual scale (z-score scale).

[Bob Carpenter]

It's actually the condition of the Hessian (ratio of largest to
smallest eigenvalue of the matrix of second derivatives of the
posterior evaluated at a particular point). What you see in
the funnel (hierarchical model with no data) is a very flat vertical
and a very sharp horizontal (at least as plotted in the manual)---that is,
a very steep ridge. That's very badly conditioned and very hard to sample
because you need increasingly small stepsizes to move the in the x direction,
which then won't let you move in the y direction.

The only recourse in these cases is to reparameterize (at least until
we get Riemannian HMC ready for prime time in the interfaces---you can
already use it in C++ directly---it's part of stan-dev/stan now).

- Bob

[George Peck]

Rethinking whether my predictors are continuous or categorical.  Trap type is categorical, and for this analysis is best left categorical because in my ms I will focus on the effect of trap type on mosquito catch rates.  Date was the time between trap deployments: 24, 48, or 72 hrs (and it should matter as an area may be "over-trapped" and need time to recover its mosquito population).  Location was a rotation among forest trap deployment sites (4 sites) and we tried to choose sites that were equal in forest cover so I would hope that inclusion of "location" in the model will have little predictive value.  I also have temperature and humidity data, and could include this in the model.  Study was done in Iquitos, Peru (deep in the Amazon rain forest), so temp and humidity fluctuated very little during study.  Thoughts?

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/kagwRgoddPY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[George Peck]

Rethinking whether my predictors are continuous or categorical.  Trap type is categorical, and for this analysis is best left categorical because in my ms I will focus on the effect of trap type on mosquito catch rates.  Date was the time between trap deployments: 24, 48, or 72 hrs (and it should matter as an area may be "over-trapped" and need time to recover its mosquito population).  Location was a rotation among forest trap deployment sites (4 sites) and we tried to choose sites that were equal in forest cover so I would hope that inclusion of "location" in the model will have little predictive value.  I also have temperature and humidity data, and could include this in the model.  Study was done in Iquitos, Peru (deep in the Amazon rain forest), so temp and humidity fluctuated very little during study.  Thoughts?

On Monday, August 22, 2016 at 9:47:53 PM UTC-7, George Peck wrote:

[Bob Carpenter]

Just make sure that if you have a group of K categorical predictors
that you only have coefficients ("intercepts") for (K - 1) of them.
Because otherwise, if you have a top-level intercept or multiple
such groups of coefficients, the model won't be identifiable.

And if there are more than a couple of group members, you can fit
a hirarchical prior.

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

[George Peck]

Have added three more predictors (rainfall, temperature, relative humidity) and realize that with six predictors the number of interactions is very large.  Is there a diagnostic to check for possibility of interactions among predictors in chapter 7 of Statistical Rethinking?  Theoretically, I expect that only the trap predictor will be important, but I would like evidence in the ms that I did test relevant predictors and their interactions.

On Wed, Aug 24, 2016 at 4:08 AM, Richard McElreath <rmcel...@gmail.com> wrote:

--

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/kagwRgoddPY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.