[R-sig-ME] wider than expected confidence intervals with lsmeans and predict.glmmadmb

[Evan Palmer-Young]

Dear List,

I am trying to use lsmeans to get confidence intervals for different levels
of treatment.

I was surprised to find that even when a fixed effect in my model was
highly significant, the confidence intervals on the lsmeans plot overlapped
almost completely. I reproduced this behavior with the "Owls" dataset. The
lsmeans() function and the predict.glmmadmb() function both gave the same
result, so there do not appear to be any surprises due to lsmeans.

I would be grateful if anybody could explain the reason for the large
confidence bands despite the significant fixed effect.


​Here is a short reproducible example-- thanks very much for any insight!

library(glmmADMB)

library(lsmeans)

#Use data from Bolker et al 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=TRUE,
         family="nbinom")
summary(m.nb)
# Estimate Std. Error z value Pr(>|z|)
# (Intercept)             4.2674     0.4705    9.07  < 2e-16 ***
#  * FoodTreatmentSatiated  -0.2602     0.0845   -3.08   0.0021 ** *
#   ArrivalTime            -0.0840     0.0190   -4.42  9.8e-06 ***
#Plot lsmeans by FoodTreatment
owls.lsm<-lsmeans(m.nb, ~FoodTreatment)
owls.lsm
# FoodTreatment   lsmean        SE df asymp.LCL asymp.UCL
# Deprived      2.188727 0.7205142 NA 0.7765454  3.600909
# Satiated      1.928499 0.7498151 NA 0.4588887  3.398110
#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.188727 0.7205142
# 2      Satiated    24.75763 1.928499 0.7498151
#Matches with lsmeans output
​


-- 
Evan Palmer-Young
PhD candidate
Department of Biology
221 Morrill Science Center
611 North Pleasant St
Amherst MA 01003
https://scholar.google.com/citations?user=VGvOypoAAAAJ&hl=en
https://sites.google.com/a/cornell.edu/evan-palmer-young/
epalmery at cns.umass.edu
ecp52 at cornell.edu

	[[alternative HTML version deleted]]