Skip to content

GLMER vs. MCMCglmm - difference in model results - Poisson model with offset term

2 messages · Bernd Lenzner, Paul Johnson

Bernd Lenzner
# 
Dear everyone,

I have a question regarding the comparison of Generalized Linear Mixed 
Models (GLMM - using the lme4 - package) and Markov-Chain Monte-Carlo 
Generalized Mixed Models (MCMCglmm - using the MCMCglmm-package).

I am running mixed models using both packages. While the estimates 
provided by each approach are relatively comparable, the significance 
levels differ strongly. My main questions are the following:

1. How can the divergence in results be explained and why doe both 
methods provide different results even though the model structure is the 
same? (Note: diagnostic plots for both models look fine)

2. I am aware that there is strong debate about the validity of p-values 
within GLMMs. Is there similar concern regarding MCMCglmm?

3. What would be a sensible interpretation of both model outputs 
especially with respect to PRED2 and the interaction term.


Below I provide a short description and example output:

I want to test the effect of 3 continuous predictors (PRED2, PRED3, 
PRED4) on a proportional response (RESP). PRED2 and PRED3 enter the 
model as an interaction effect (PRED2*PRED3). The response is modeled 
using a Poisson error structure with a log-link function and the 
nominator of the ratio is therefore included as an offset predictor 
(offset(log(PRED1))).


MODEL1 -  GLMM using the lme4-package:


MOD1 <- glmer(RESP ~ offset(log(PRED1)) + PRED2 * PRED3 + PRED4 + 
(1|RANEF), family = "poisson", data = data, 
control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000)))

summary(MOD1)
Generalized linear mixed model fit by maximum likelihood (Laplace 
Approximation) ['glmerMod']
  Family: poisson  ( log )
Formula: RESP ~ offset(log(PRED1)) + PRED2 * PRED3 + PRED4 + (1 | RANEF)
    Data: data
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

      AIC      BIC   logLik deviance df.resid
   2451.3   2470.9  -1219.6   2439.3      190

Scaled residuals:
      Min       1Q   Median       3Q      Max
-10.8175  -1.2477  -0.1913   1.1076  11.0389

Random effects:
  Groups Name        Variance Std.Dev.
  RANEF  (Intercept) 0.1519   0.3898
Number of obs: 196, groups:  Order, 56

Fixed effects:
                   Estimate    Std. Error z value Pr(>|z|)
(Intercept)       -3.29896    0.06250    -52.784   < 2e-16 ***
PRED2              0.09089    0.01850     4.912     9.00e-07 ***
PRED3              0.57262    0.02375    24.114    < 2e-16 ***
PRED4              0.54726    0.01808    30.275    < 2e-16 ***
PRED2:PRED3        0.12190    0.01821     6.693     2.19e-11 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

Correlation of Fixed Effects:
               (Intr)    PRED2     PRED3    PRED4
PRED2         -0.142
PRED3         -0.056    0.102
PRED4         -0.042    0.135    -0.155
PRED2:PRED3    0.001    0.017    -0.678     0.273


MODEL2 -  MCMCglmm using the lme4-package:

Now I build the same model using a MCMCglmm using weak informative 
priors. The offset term was implemented following the suggestion by 
Jarrod Hadfield in this post 
(https://stat.ethz.ch/pipermail/r-sig-mixed-models/2017q3/026016.html). 
See the model formula and summary below:


prior <- list(B=list(V=diag(6)*10, mu=c(0,1,0,0,0,0)),
                          G =list(G1=list(V=1,nu=0.002)),
                          R=list(V=1,nu=0.002))
prior$B$V[2,2] <- 1e-9 # set fixed effect prior for the offset variable to 1


PGLMM.MOD1 <- MCMCglmm(RESP ~ log(PRED1) + PRED2 * PRED3 + PRED4,
                                         random = ~ RANEF,
                                         family="poisson",
                                         prior = prior,
                                         scale = T,
                                         pr=TRUE,
                                         data = data,
                                         nitt = 500000, burnin = 20000, 
thin = 200)


summary(PGLMM.MOD1)

Iterations = 20001:499801
  Thinning interval  = 200
  Sample size  = 2400

  DIC: 1164.591

  G-structure:  ~Order

         post.mean  l-95% CI    u-95% CI   eff.samp
RANEF   0.02264    0.0002387   0.0762      2400

  R-structure:  ~units

          post.mean    l-95% CI     u-95% CI     eff.samp
units    0.2788       0.1946       0.3705        2400

  Location effects: RESP ~ log(PRED1) + PRED2 * PRED3 + PRED4

              post.mean   l-95% CI  u-95% CI   eff.samp   pMCMC
(Intercept)  -3.40280   -3.50958   -3.28826   2400       <4e-04 ***
log(PRED1)   1.00000    0.99994    1.00006    2257       <4e-04 ***
PRED2        0.08369    -0.01995   0.18120    2400       0.113
PRED3        0.63058    0.52164    0.73889    2400       <4e-04 ***
PRED4        0.60488    0.50565    0.69963    2400       <4e-04 ***
PRED2:PRED3  0.06882    -0.02363   0.17414    2258       0.167
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1


As you can see from the different summaries, the estimates are fairly 
similar with some divergence especially regarding PRED2 and the 
interaction term. Additionally the provided p-values differ 
dramatically, with the GLMM (lme4) being highly significant for PRED2 
and the interaction wheras non-significance for the MCMCglmm.

My questions now are:

1. How can the divergence in results be explained and why doe both 
methods provide different results even though the model structure is the 
same? (Note: diagnostic plots for both models look fine)

2. I am aware that there is strong debate about the validity of p-values 
within GLMMs. Is there similar concern regarding MCMCglmm?

3. What would be a sensible interpretation of both model outputs 
especially with respect to PRED2 and the interaction term.


Any insight and help on that matter would be highly appreciated.

Best regards,

Bernd
4 days later
Paul Johnson
# 
Hi Bernd,

I think there are two problems. 

1. 
If your response is a number of successes (RESP) out of a number of trials (PRED1), it should be modelled with the binomial distribution. In glmer this would take the form cbind(RESP, PRED1 - RESP) ~ ? . (If the proportion is low, you?ll likely find this makes little difference.)

2. 
The two Poisson models aren?t comparable because your glmer Poisson model has no residual term, i.e. no random effect to mop up variation not explained by the fixed effects and the Poisson distribution, while MCMCglmm fits this by default. You can see this in the glmer output in the residual deviance being much larger than the residual df, and in the wide spread of the standardised residuals. The basic problem is that a pure Poisson GLMM assumes that all the variation in your model can be explained, except for the very limited degree of random variation allowed by the Poisson distribution. This is hardly ever a sensible starting assumption when modelling messy biological count data (which is what most of us are doing). It?s a very common error and quite dangerous because it tends to generate spuriously significant results, as you appear to have found. The same applies to comparing binomial models (provided n trials > 1) between glmer and MCMCglmm. If you add an observation-level random intercept to the glmer model (whether binomial or Poisson) then the two models will be comparable, i.e. (1 | obs) where obs is factor(1:nrow(data)). You can also model overdispersion with a negative binomial using glmer.nb, but there?s no comparable method in MCMCglmm.

Al the best,
Paul 

PS as an aside, you can spot overdispersion in the standard residuals-by-fitted diagnostic scatter plot by comparing it with a plot simulated from the fitted model:

library(devtools)
install_github("pcdjohnson/GLMMmisc")
library(GLMMmisc)
library(lme4)

# Poisson-lognormal model with random effect to model overdispersion (INDEX)
fit.poisln <-
  glmer(TICKS ~ YEAR + scale(HEIGHT) + (1 | BROOD) + (1 | LOCATION) + (1 | INDEX),
        family = "poisson", data = grouseticks)
# pure Poisson model with no random effect to model overdispersion
fit.pois <-
  glmer(TICKS ~ YEAR + scale(HEIGHT) + (1 | BROOD) + (1 | LOCATION),
        family = "poisson", data = grouseticks)

par(mfrow = c(1, 2))
sim.residplot(fit.pois)
title("Pure Poisson GLMM")
sim.residplot(fit.poisln)
title("Poisson-lognormal GLMM")

# if you run the plotting lines a few times you'll see that the simulated residuals 
# tend to look simular to the real residuals from the Poisson-lognormal GLMM
# (a good sign) but much less spread out than the real residuals from the 
# pure Poisson model, a sign that the latter model is failing to account for
# a substantial amount of variation in the responses.
# NB the DHARMa package does something similar but in a more formal and
# sophisticated way