[R-sig-ME] GLMER vs. MCMCglmm - difference in model results - Poisson model with offset term

Bernd Lenzner bernd@|enzner @end|ng |rom un|v|e@@c@@t
Fri May 24 10:51:26 CEST 2019


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


-- 
Bernd Lenzner PhD

Researcher
University of Vienna
Department of Botany and Biodiversity Research
Division of Conservation Biology, Vegetation and Landscape Ecology
Rennweg 14
A-1030 Vienna



More information about the R-sig-mixed-models mailing list