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

Martin Maechler m@ech|er @end|ng |rom @t@t@m@th@ethz@ch
Fri May 24 09:58:40 CEST 2019


Dear Bernd,

can you please repost using "plain text" instead of "rich text"
/ "formatted text" or whatever your e-mail app calls what really
is HTML ?

See how almost unreadable your  R output  ends up in the post of
your e-mail :

>>>>> Bernd Lenzner 
>>>>>     on Thu, 23 May 2019 16:50:23 +0200 writes:

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


    > 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


    > [[alternative HTML version deleted]]

    > _______________________________________________
    > R-sig-mixed-models using r-project.org mailing list
    > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models



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