[R-sig-ME] wider than expected confidence intervals with lsmeans and predict.glmmadmb
John Maindonald
john.maindonald at anu.edu.au
Sat May 27 03:02:23 CEST 2017
The models m.nb2 (from fitting using glmmTMB()) and m.nb (from glmmadmb())
return coefficient and SE information that is for all practical purposes identical
> m.nb2$call
glmmTMB(formula = NCalls ~ FoodTreatment + ArrivalTime + +(1 |
Nest), data = Owls, family = "nbinom2", ziformula = ~0, dispformula = ~1)
> print(coef(summary(m.nb2)), digits=2)
$cond
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.91 0.633 7.8 9.1e-15
FoodTreatmentSatiated -0.69 0.107 -6.5 9.4e-11
ArrivalTime -0.12 0.025 -4.6 4.9e-06
> m.nb$call
glmmadmb(formula = NCalls ~ FoodTreatment + ArrivalTime + +(1 |
Nest), data = Owls, family = "nbinom", zeroInflation = FALSE)
> print(coef(summary(m.nb)), digits=2)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.91 0.633 7.8 9.1e-15
FoodTreatmentSatiated -0.69 0.107 -6.5 9.4e-11
ArrivalTime -0.12 0.025 -4.6 4.9e-06
The differences between the two graphs are then a conseqence of
what is done on the way to creating those graphs. For comparing
levels of `FoodTreatment`, the SE is for the in each case for the
difference, not for the levels individually. It is then information that it
would be helpful to add to the graph given by:
owls.lsm<-lsmeans(m.nb, ~FoodTreatment)
plot(owls.lsm)
One can get a plot that shows the SE for the difference thus:
K <- diag(length(coef(m.nb)))[1:2,]
rownames(K) <- c("Deprived","Sat-Dep”)
library(multcomp)
plot(glht(m.nb,linfct=K))
Or, nearer to what you want, maybe:
K2 <- rbind(K[1,], c(1,1,0), K[2,])
rownames(K2) <- c("Deprived","Saturated", "Sat-Dep")
plot(glht(m.nb,linfct=K2))
It is, of course, in this simple case, possible to place intervals
around the two estimates, designed so that if the intervals do
not overlap, then the difference is not “significant” at alpha=0.05.
I will leave it to you, or to others, to check just what the code you
give, that uses as its starting-point output from glmmTMB(), may be
doing. This is not a straightforward use of lsmeans().
John Maindonald email: john.maindonald at anu.edu.a<mailto:john.maindonald at anu.edu.a>
On 27/05/2017, at 09:29, Evan Palmer-Young <ecp52 at cornell.edu<mailto:ecp52 at cornell.edu>> wrote:
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<mailto: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<mailto:john.maindonald at anu.edu.au>.
On 26/05/2017, at 01:04, Evan Palmer-Young <ecp52 at cornell.edu<mailto: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<mailto:epalmery at cns.umass.edu>
ecp52 at cornell.edu<mailto:ecp52 at cornell.edu>
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-models at r-project.org<mailto: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<mailto:epalmery at cns.umass.edu>
ecp52 at cornell.edu<mailto:ecp52 at cornell.edu>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list