[R-sig-ME] Repost in plain text: Discrepancy in residuals in equivalent glmer (lme4) and glmmadmb models

[Adam Smith]

Sorry all, here's a plain text version that is hopefully more user-friendly...

Original post follows:

All,

Can someone please clarify for me why lme4 (glmer) and glmmADMB would generate essentially identical estimate of fixed and random effects in a Poisson log-normal (individual-level random effect for overdispersion) but give (often dramatically) different residual values?  It must have something to do with the way fitted values are estimated or SDs are calculated for each observation, yes? 

Any reason I should favor one over the other?

For context, I was comparing the Poisson log-normal model with negative binomial models in glmmADMB, but I also like to compare the lme4 estimation with the glmmADMB version for peace of mind.

Thanks for any clarifications...

Adam Smith

# Loading required packages
toLoad <- c("glmmADMB", "lme4", "ggplot2", "gridExtra")
lapply(toLoad, library, character.only = TRUE)

> sessionInfo()
R version 2.14.1 (2011-12-22)
Platform: x86_64-pc-mingw32/x64 (64-bit)

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets  methods   base     

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] gridExtra_0.9.1    ggplot2_0.9.1      numDeriv_2010.11-1 lme4_0.999375-42   Matrix_1.0-3       lattice_0.20-0    
[7] glmmADMB_0.7.2.5   R2admb_0.7.5       MASS_7.3-16       

# Data available here if you feel like messing with it:
# https://www.dropbox.com/s/9zzy5qb8y207os0/testdata
# NOTE: When downloaded, I think dropbox will attach a*.txt extension to the file
# Delete the extension before loading with dget()
testdata = dget("C:/temp/testdata") # Obviously, code assume you saved it to C:/temp
str(testdata)

# Code for lme4 and glmmADMB versions of the Poisson log-normal GLMM
oP_glmer <- glmer(passes ~ phen + nightWsp + nightMb + nightWp12 + nightWpSeas + deltaTemp +
  deltaRh + deltaMb + propPrecip + tempresid + skycov + vis + (1|site:year) + (1|obs), 
                  family="poisson", data=testdata)

oP_admb <- glmmadmb(passes ~ phen + nightWsp + nightMb + nightWp12 + nightWpSeas + deltaTemp +
  deltaRh + deltaMb + propPrecip + tempresid + skycov + vis + (1|site:year) + (1|obs),
                    family="poisson", data=testdata, admb.opts=admbControl(shess=FALSE,noinit=FALSE))

# Creating plots comparing fixed and random effects between models
fixed.df <- data.frame(glmer_fixed = fixef(oP_glmer), admb_fixed = fixef(oP_admb))
fixed_plot <- ggplot(fixed.df, aes(glmer_fixed, admb_fixed)) + geom_smooth() + geom_point()
site_year.df <- data.frame(ranef(oP_glmer)[[2]][,1], ranef(oP_admb)[[1]]); names(site_year.df) = c("glmer_re1", "admb_re1")
site_year_plot <- ggplot(site_year.df, aes(glmer_re1, admb_re1)) + geom_smooth() + geom_point()
obs.df <- data.frame(ranef(oP_glmer)[[1]][,1], ranef(oP_admb)[[2]]); names(obs.df) = c("glmer_re2", "admb_re2")
obs_plot <- ggplot(obs.df, aes(glmer_re2, admb_re2)) + geom_smooth() + geom_point()

# Creating plot comparing residuals between models
fitted.df <- data.frame(glmer_fitted = fitted(oP_glmer), admb_fitted = fitted(oP_admb))
fitted_plot <- ggplot(fitted.df, aes(glmer_fitted, admb_fitted)) + geom_smooth() + geom_point()
resids.df <- data.frame(glmer_resids = residuals(oP_glmer, type="pearson"), admb_resids = residuals(oP_admb, type = "pearson"))
resid_plot <- ggplot(resids.df, aes(glmer_resids, admb_resids)) + geom_smooth() + geom_point()

# Generating plot
grid.arrange(arrangeGrob(fixed_plot, site_year_plot, obs_plot, ncol = 3), fitted_plot, resid_plot)