[R-sig-ME] GLMER vs. MCMCglmm - difference in model results - Poisson model with offset term
Paul Johnson
p@u|@john@on @end|ng |rom g|@@gow@@c@uk
Tue May 28 10:58:25 CEST 2019
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
> On 24 May 2019, at 09:51, Bernd Lenzner <bernd.lenzner using univie.ac.at> wrote:
>
>
> 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
>
> _______________________________________________
> 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