help with Gamma GLMM inconsistent results between nimble and glmer implementations

[Tiago Marques]

Dear NIMBLE users group,

 

Sorry upfront for the long question, please bear with me. I have implemented in Nimble a Bayesian equivalent to a Gamma GLMM model (see model code below) previously fit in lme4:glmer

 

glmer(rate~(1|location)+(1|year),data=data,family=Gamma(link="log"))

 

where both factors, location and year are random effects. Even though mixing and convergence seem ok I am having trouble reproducing the results from glmer, , in particular the location random effect standard deviation and even the intercept. I know that Bayes and frequentist implementations do not need to produce the exact same results, but I used what I hoped were priors vague enough such that they should match well. Not being very experienced with Bayesian implementations these seem off enough to make me think there might be something that I am missing.

 

Could perhaps someone more experienced in Nimble (or Bayesian implementations in general) comment if these differences are indeed somewhat suspicious, and if so, comment on what might be the issue? In particular, I wonder if (cf. model code in PS1 below)

1. there is perhaps something wrong with the Gamma parameterization?

2. there is something different I should try for priors on the random effects sd/variance?

Posterior plots with glmer on top for comparison here:

https://github.com/TiagoAMarques/Report4BB/blob/main/LookingAtNimble/PosteriorDistributionsGammaGLMM.jpeg

 

Any pointers that might help resolve this question would be greatly appreciated and might in fact be repaid in the drinks of your choosing during your next visit to Lisbon.

 

best,

 

Tiago

 

PS1 - The Nimble model code looks like this:

 

  ## define the model

  GLMMcode <- nimbleCode({

    # the overall intercept

    beta0 ~ dnorm(0,sd=10)

    # random effect standard deviation associated with location, a uniform, might change this to be something else latter

    sigmal_RE ~ dunif(0, 2)

    # random effect standard deviation associated with year, a uniform, might change this to be something else latter

    sigmay_RE ~ dunif(0, 2)

    # the gamma dispersion (or variance - see commented parametrization 1) parameter, a uniform, might change this to be something else latter

    #dispersion ~ dunif(0, 10)

    disp ~ dunif(0, 10)

    #get year random effects

    for(yy in 1:nyears){

      #REy[yy] ~ dnorm(0, sd = sigmay_RE)

      REy[yy] ~ dnorm(0, sd = 1)

    } 

    #get location random effects

    for(ll in 1:nlocs){

      #REl[ll] ~ dnorm(0, sd = sigmal_RE)

      REl[ll] ~ dnorm(0, sd = 1)

    } 

    for (i in 1:N){

      #get the linear predictor, consider a log link function

      #log(mean[i]) <- beta0 + REy[year[i]] + REl[loc[i]]

      # now Using decentered parametrization, a suggestion by Ben Augustine

      log(mean[i]) <- beta0 + REy[year[i]]*sigmay_RE + REl[loc[i]]*sigmal_RE

      #parametrization 1

      # crate[i] ~ dgamma(shape=(mean[i]^2)/disp,scale=disp/mean[i])

      #parametrization 2

      crate[i] ~ dgamma(shape=1/disp,scale=mean[i]*disp)

    }

  })

Trace plots here: https://github.com/TiagoAMarques/Report4BB/blob/main/LookingAtNimble/TracePlotsGammaGLMM.jpeg

A full working data + code package is here

https://github.com/TiagoAMarques/Report4BB/tree/main/LookingAtNimble

and runs in under a minute if you are feeling especially keen to help (thank you in advance if that is the case!): just source in R the file “GLMMinNimble.R”

 

PS2 - one can never say thanks enough to those who help you. I am not a Bayesian and even the simple code above actually already results form kind and useful interactions with both Ben Augustine and Perry de Valpine. This question is coming via Ben Bolker having set me in the Bayesian way via an answer on stack exchange about how to estimate the variance on predictions for unobserved levels of a Gamma GLMM here (https://stats.stackexchange.com/questions/616697/how-to-estimate-precision-on-a-prediction-for-a-glmm-for-an-unobserved-level-of/).

[Wei Zhang]

Hi Tiago,

I played with your model a bit today and it seems to me that the Gamma parameterization in nimble does not have problems. I think one reason for the difference between nimble MCMC and lme4 is that they do not parametrize the model in exactly the same way. I cannot find how lme4 does this exactly, but I tried using glmmTMB to fit the nimble model and can confirm that glmmTMB and the nimble model you got are doing the same thing. Comparing the MLEs from glmmTMB and nimble Laplace (they are pretty close) to MCMC, the gaps seem to be smaller, although may not be as close as you hope to see. I added some code to yours and it is attached here. You can have a  quick run to see what you think. I guess using Laplace might be a reason for the difference, which anyway involves approximations. I am very happy to discuss further. 

Best wishes,

Wei

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAFOBQAq-X5_1WUAUB4Bn2%3DgFCZp1k5z61%3D8ZJODyjzFNTWGSFg%40mail.gmail.com.