residfitted.olre: Output "corrected" Pearson residuals and fitted values

In pcdjohnson/GLMMmisc: Miscellaneous functions for GLMMs

Description

Output "corrected" Pearson residuals and fitted values from a Poisson GLMM with an observation-level random effect (OLRE) fitted using lme4::glmer. Defaults to residuals(., type = "pearson") and fitted(.) for other classes of fit.

Usage

residfitted.olre(object, warn = TRUE)

Arguments

object	Fitted model object
warn	Warn when not Poisson glmer object with OLRE

Examples

# fit three models to the grouse ticks count data (see ?grouseticks)
# first, a Poisson GLMM that doesn't account for overdispersion
form0 <- TICKS ~ YEAR + scale(HEIGHT) + (1 | BROOD) + (1 | LOCATION)
fit.pois  <-
  lme4::glmer(form0, family = "poisson", data = lme4::grouseticks)

# second, model overdispersion by adding an OLRE: "+ (1 | INDEX)"
# giving a Poisson lognormal GLMM (see http://www.ncbi.nlm.nih.gov/pubmed/11393830)
# also see Xavier Harrison's article on this: https://peerj.com/articles/616/
form1 <- TICKS ~ YEAR + scale(HEIGHT) + (1 | BROOD) + (1 | LOCATION) + (1 | INDEX)
fit.poisln <-
  lme4::glmer(form1, family = "poisson", data = lme4::grouseticks)

# finally, model overdispersion with a negative binomial distribution
#fit.nb <- glmmadmb(form0, family = "nbinom", data = lme4::grouseticks)
# I'VE COMMENTED OUT LINES REQUIRING glmmADMB AS THEY WERE CAUSING PROBLEMS
# WITH BUILDING THE PACKAGE -- I'LL TRY TO FIX THIS

# calculate Pearson residuals and fitted values
residfitted.olre(fit.pois)
#residfitted.olre(fit.nb)
residfitted.olre(fit.poisln)
# if the model isn't a Poisson-lognormal GLMM the standard stats functions are used
# with an optional warning.

# both the Poisson-lognormal and the NB fit much better than the Poisson,
# and the NB fits slightly better than the Poisson-lognormal
#AIC(fit.pois, fit.poisln, fit.nb)

# assess the fit of the model by plotting Pearson residuals against fitted values
par(mfrow = c(2, 2))
plotloess(x = fitted(fit.pois), y = residuals(fit.pois, type = "pearson"))
title("Poisson GLMM residuals v fitted plot")
plotloess(x = fitted(fit.poisln), y = residuals(fit.poisln, type = "pearson"))
title("Poisson-lognormal GLMM residuals v fitted plot")
#plotloess(x = fitted(fit.nb), y = residuals(fit.nb, type = "pearson"))
#title("Negative binomial GLMM residuals v fitted plot")

# the standard residuals v fitted plots look fine for the Poisson & NB fits,
# but the Poisson-lognormal plot shows a nasty trend and severe heteroscedasticity.
# this is because the overdispersion effects are included in the fitted values.
# to be comparable with the negative binomial GLMM plot, the residuals need to be
# subtracted from the fitted values and added to the residuals. this is what the
# residfitted.olre function does:
plotloess(residfitted.olre(fit.poisln))
title("Corrected Poisson-lognormal GLMM residuals v fitted plot")

# now the residuals v fitted plot from the Poisson-lognormal GLMM can be
# assessed alongside the other two models using the usual criteria: linearity
# and homoscedasticity. the residuals from the NB and Poisson-lognormal distributions
# look very similar, which isn't surprising because each is a Poisson distribution
# with added (really mutliplied) gamma- and lognormal-distributed noise
# respectively, and gamma and lognormal distributions can be quite similar.

pcdjohnson/GLMMmisc documentation built on May 24, 2019, 11:40 p.m.