[R-sig-ME] GoodnessOfFit-and-PseudoR²computingFor-glmmTMB-models

[Paul Johnson]

No problem. I didn’t read your questions properly though, apologies if part of my reply didn’t make sense. I'll try to answer questions 1 and 3 as well:

Question 1. 
A good way is to plot your data over your predictions and visually assess fit. This isn’t always easy, especially if you have multiple continuous predictors, in which case you can produce multiple predicted lines conditioned on a few values. Another way is to plot the Pearson residuals against fitted values, and check for homoscedasticity along the fitted values axis, and absence of any nonlinear trend. A good way is to simulate response data from the fitted model (which must fit perfectly) a few times, and compare it with the real data using the methods above (plotting data vs predictions, residuals vs fitted values plot). If the real data fits the model well, it should look similar to the simulated data.

Question 3.
I’ve pasted below a script showing how to estimate R2 for a Poisson GLMM adapting one of the examples from glmmTMB, following this paper: 
"The coefficient of determination R2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded"
http://rsif.royalsocietypublishing.org/content/14/134/20170213
which built on 
“A general and simple method for obtaining R2 from generalized linear mixed‐effects models”
doi.org/10.1111/j.2041-210x.2012.00261.x

There are a couple of packages that have functions to do this (piecewiseSEM, MuMIn), but they don’t work on glmmTMB (yet — piecewiseSEM is actively maintained so this might change).

All the best,
Paul

# add observation-level factor to Salamanders data
Salamanders$obs <- factor(1:nrow(Salamanders))
(m1 <- glmmTMB(count~ mined + (1|site) + (1|obs), 
               family=poisson, data=Salamanders))
summary(m1)

# fixed effects variance (manually - not sure how to get 
# predict.glmmTMB to do this)
linear.predictor <- model.matrix(m1) %*% fixef(m1)$cond
fixed.var <- var(linear.predictor)
# sum variance of all random effects ***excluding OLRE***
all.ranef.var <- unlist(VarCorr(m1)$cond)
ranef.var <- all.ranef.var[!names(all.ranef.var) %in% "obs"]
# OLRE (additive overdispersion variance)
olre.var <- all.ranef.var["obs"]

# now the observation-level variance 
# (distinct from the observation-level random effect [OLRE])
# here this is the variance added by the Poisson distribution 
# (aka distribution-specific variance)

# first we need to estimate lambda
# the Interface paper recommends calculating lambda as 
# exp(beta0 + total.re.var/2)   -- eqn. 5.8
# beta0 is the intercept from an intercept-only refit of the model
(beta0 <- fixef(update(m1, ~ . - mined))$cond)
lambda <- exp(beta0 + (ranef.var + olre.var)/2)

# observation-level variance (Table 1)
ol.var <- trigamma(lambda) 

# total variance
total.var <- fixed.var + ranef.var + olre.var + ol.var

# marginal R2glmm
fixed.var / total.var
# conditional R2glmm
(fixed.var + ranef.var) / total.var



> On 6 Feb 2018, at 22:13, C. AMAL D. GLELE <altessedac2 at gmail.com> wrote:
> 
> Many thanks, Paul.
> Regards,
> 
> 2018-02-06 15:53 GMT+01:00 Paul Johnson <paul.johnson at glasgow.ac.uk>:
> Hi,
> 
> My preferred approach to overdispersion is none of 1-3 but to assume it applies (it usually does, for biological data anyway), and make sure the model includes a parameter to model overdispersion. For binomial (except Bernoulli) and Poisson you can include an observation-level random effect (OLRE). You can then gauge the amount of overdispersion in the model from the size of the OLRE variance estimate. (Note the OLRE will mop up variation due to *all* sources of lack of fit, including poor model specification, e.g. fitting a straight line where a curve would fit better.)
> 
> Best wishes,
> Paul
> 
> 
> > On 6 Feb 2018, at 14:29, C. AMAL D. GLELE <altessedac2 at gmail.com> wrote:
> >
> > Hi, dear all
> > Please, your help for the following problem:
> > when I fit a mixed model using glmmTMB (poisson family or others), how do I:
> > 1) check, if the model fits well my data (goodness of fit)?
> > 2) check if my model is overdisperced or not (by using sigma(model)?)
> > 3) compute an pseudo R² to see the percentage of the variability of my
> > response which is explained by my model?
> > In advance, thanks for your answers.
> > Regards,
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
>