[R-sig-ME] wider than expected confidence intervals with lsmeans and predict.glmmadmb
Evan Palmer-Young
ecp52 at cornell.edu
Fri May 26 23:29:52 CEST 2017
Thanks very much for your reply, Prof. Maindonald.
I agree that the pairwise comparisons are informative, but it would be
easiest for readers to see the data on the original scale to show
differences between groups.
When the lsmeans are plotted from glmmTMB, which fits a model with fixed
effects identical to those in glmmADMB, the estimates are identical but the
SE's differ by a factor of 8.
So I am still confused about why the lsmeans plots would reflect pairwise
differences with some packages but not with glmmADMB.
In my experience, lsmeans plots of group means from glmer() models are also
non-overlapping when pairwise comparisons are highly significant.
I have extended the code to illustrate the differences.
library(glmmADMB)
library(lsmeans)
#Use data from worked example
#http://glmmadmb.r-forge.r-project.org/glmmADMB.html
library(glmmADMB)
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.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
################## Compare 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 #######
#Ben Bolker wrote a function to talk to lsmeans-- incredible!
# 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 ###
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?!
Thank you again.
Evan
On Thu, May 25, 2017 at 5:22 PM, John Maindonald <john.maindonald at anu.edu.au
> wrote:
> The confidence intervals that you have obtained are for the levels
> of `FoodTreatment`, not for the contrast `Satiated-Deprived`.
>
>
> Try, the following, which also gives a confidence interval for the
> difference from the initial level of `FoodTreatment`:
>
> > library(glmmADMB)
> > library(lsmeans)
> > 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”)
> > owls.lsm<-lsmeans(m.nb, ~FoodTreatment)
> > lsmeans (owls.lsm, "FoodTreatment", contr = "trt.vs.ctrl")
> $lsmeans
> . . .
>
> $contrasts
> contrast estimate SE df z.ratio p.value
> Satiated - Deprived -0.260228 0.084501 NA -3.08 0.0021
>
> John Maindonald email: john.maindonald at anu.edu.au.
>
>
> On 26/05/2017, at 01:04, Evan Palmer-Young <ecp52 at cornell.edu> wrote:
>
> 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]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
>
--
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]]
More information about the R-sig-mixed-models
mailing list