predictions and lsmeans: Very large standard errors compared to other models with same fixed effects

[ecp52]

I have a model fit in glmmADMB for which a predictor variable has a significant effect. When I plot prediction intervals for equivalent models in glmmTMB and lme4, the prediction intervals are consistent with a strong fixed effect, as indicated by the model summary. However, the prediction intervals derived from glmmADMB are about 8 times as wide, despite the fact that the fixed effects coefficients are identical in glmmADMB and glmmTMB models, and pariwise comparisons between groups are also the same.

I cannot figure out why these discrepancies are so large. I post this to the glmmadmb page because the standard errors for glmmadmb seem to be the anomalies-- those from glmmTMB and lme4 were similar.
The glmmTMB and lme4 prediction intervals also seem more plausible than the glmmADMB intervals, i.e., there are non-overlapping confidence intervals for significantly different groups.

Any suggestions or explanations would be greatly appreciated. I've been very grateful to use glmmADMB for its negative binomial distribution and crossed random effects, and would like to create plots that accurately reflect the model predictions. Thanks!

Here is code to reproduce the phenomenon:

#install.packages("glmmADMB", repos = "http://glmmadmb.r-forge.r-project.org/repos")

library(glmmADMB)
library(glmmTMB)
library(lsmeans)

#Use data from worked example
#http://glmmadmb.r-forge.r-project.org/glmmADMB.html

data(Owls)
str(Owls)
Owls <- transform(Owls, 
                  Nest=reorder(Nest,NegPerChick), 
                  logBroodSize=log(BroodSize), 
                  NCalls=SiblingNegotiation)


m.nb<- glmmadmb(NCalls~FoodTreatment+ArrivalTime+ 
           +(1|Nest), 
         data=Owls, 
         zeroInflation=FALSE, 
         family="nbinom")
summary(m.nb)
# Estimate Std. Error z value Pr(>|z|)    
# (Intercept)             4.9101     0.6334    7.75  9.1e-15 ***
#   FoodTreatmentSatiated  -0.6924     0.1069   -6.48  9.4e-11 ***
#   ArrivalTime            -0.1154     0.0253   -4.57  4.9e-06 ***
#Plot lsmeans by FoodTreatment 
owls.lsm<-lsmeans(m.nb, ~FoodTreatment)
owls.lsm 
# FoodTreatment   lsmean        SE df  asymp.LCL asymp.UCL
# Deprived      2.053073 0.8952071 NA  0.2984988  3.807646
# Satiated      1.360690 0.9037320 NA -0.4105918  3.131973
#SE is much higher than for fixed effects in model

plot(owls.lsm)
 #95% confidence bands overlap almost entirely

#Confirm with predict.glmmadmb:
New.data<-expand.grid(FoodTreatment= levels(Owls$FoodTreatment),
                      ArrivalTime = mean(Owls$ArrivalTime))

New.data$NCalls <- predict(m.nb, New.data, re.form=NA, SE.fit = TRUE)

#Get standard errors:
calls.pred<- predict(m.nb, New.data, re.form = NA, se.fit = TRUE)
calls.pred<-data.frame(calls.pred)

New.data$SE<-calls.pred$se.fit
New.data
# FoodTreatment ArrivalTime   NCalls        SE
# 1      Deprived    24.75763 2.053073 0.8952071
# 2      Satiated    24.75763 1.360690 0.9037320
#Matches with lsmeans output



##################  Compare glmmADMB fit to glmmTDMB  ####################
#install.packages("glmmTMB")
library(glmmTMB)
m.nb2<- glmmTMB(NCalls~FoodTreatment+ArrivalTime+ 
                  +(1|Nest), 
                data=Owls, 
                family="nbinom2")
summary(m.nb2)
# Estimate Std. Error z value Pr(>|z|)    
# (Intercept)            4.91011    0.63343   7.752 9.07e-15 ***
#   FoodTreatmentSatiated -0.69238    0.10692  -6.476 9.44e-11 ***
#   ArrivalTime           -0.11540    0.02526  -4.569 4.90e-06 ***

#Compare to glmmADMB model:Fixed effects are identical
summary(m.nb)
# Estimate Std. Error z value Pr(>|z|)    
# (Intercept)             4.9101     0.6334    7.75  9.1e-15 ***
#   FoodTreatmentSatiated  -0.6924     0.1069   -6.48  9.4e-11 ***
#   ArrivalTime            -0.1154     0.0253   -4.57  4.9e-06 ***

#Plot lsmeans by FoodTreatment 
owls.lsm<-lsmeans(m.nb2, ~FoodTreatment)
#oops, lsmeans can't use glmmTMB object!

  ########   Interlude   #######
#Use Ben Bolker's function to talk to lsmeans
# https://github.com/glmmTMB/glmmTMB/issues/205
recover.data.glmmTMB <- function(object, ...) {
  fcall <- getCall(object)
  recover.data(fcall,delete.response(terms(object)),
               attr(model.frame(object),"na.action"), ...)
}
lsm.basis.glmmTMB <- function (object, trms, xlev, grid, vcov.,
                               mode = "asymptotic", component="cond", ...) {
  if (mode != "asymptotic") stop("only asymptotic mode is available")
  if (component != "cond") stop("only tested for conditional component")
  if (missing(vcov.)) 
    V <- as.matrix(vcov(object)[[component]])
  else V <- as.matrix(.my.vcov(object, vcov.))
  dfargs = misc = list()
  if (mode == "asymptotic") {
    dffun = function(k, dfargs) NA
  }
  ## use this? misc = .std.link.labels(family(object), misc)
  contrasts = attr(model.matrix(object), "contrasts")
  m = model.frame(trms, grid, na.action = na.pass, xlev = xlev)
  X = model.matrix(trms, m, contrasts.arg = contrasts)
  bhat = fixef(object)[[component]]
  if (length(bhat) < ncol(X)) {
    kept = match(names(bhat), dimnames(X)[[2]])
    bhat = NA * X[1, ]
    bhat[kept] = fixef(object)[[component]]
    modmat = model.matrix(trms, model.frame(object), contrasts.arg = contrasts)
    nbasis = estimability::nonest.basis(modmat)
  }
  else nbasis = estimability::all.estble
  list(X = X, bhat = bhat, nbasis = nbasis, V = V, dffun = dffun, 
       dfargs = dfargs, misc = misc)
}

#####   End interlude ###

#### Extract lsmeans from glmmTMB with the helper function ###

lsm.TMB<- lsmeans(m.nb2, ~FoodTreatment)
plot(lsm.TMB)  #non-overlapping CI's

#Compare SE's
owls.lsm
# FoodTreatment   lsmean        SE df  asymp.LCL asymp.UCL
# Deprived      2.053073 0.8952071 NA  0.2984988  3.807646
# Satiated      1.360690 0.9037320 NA -0.4105918  3.131973

lsm.TMB
# FoodTreatment   lsmean        SE df asymp.LCL asymp.UCL
# Deprived      2.053065 0.1068562 NA  1.843631  2.262500
# Satiated      1.360683 0.1161322 NA  1.133068  1.588298

#lsmeans are identical but SE's differ by factor of 8?!

#Confirm with predict()
New.data2<-expand.grid(FoodTreatment= levels(Owls$FoodTreatment),
                      ArrivalTime = mean(Owls$ArrivalTime) )

calls.pred2 <- predict(m.nb2, New.data2, re.form=NA, SE.fit = TRUE)
# Error in predict.glmmTMB(m.nb2, New.data2, re.form = NA, SE.fit = TRUE) : 
#   re.form not yet implemented


 
##Compare to lme4 fit:
library(lme4)

nb4<-glmer.nb(NCalls~FoodTreatment+ArrivalTime+ 
                +(1|Nest), 
              data=Owls)

#Convergence warning
nb4.lsm<-lsmeans(nb4, ~FoodTreatment)

plot(nb4.lsm) #well-separated, glmmADMB SE's seem to be anomalous?



######### Follow-up: ##########################

##Check Russell Lenth's function to extract data from glmmADMB
#http://permalink.gmane.org/gmane.comp.lang.r.lme4.devel/12053
#Similar to Ben Bolker's function to extract from glmmTMB
#https://github.com/glmmTMB/glmmTMB/issues/205

#I cannot figure out why they give comparable results for pairwise comparisons,
  #but different results for standard errors

[bbolker]

For reference, quoting discussion from John Maindonald on r-sig-mixed-models: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2017q2/025698.html