wider than expected confidence intervals with lsmeans and predict.glmmadmb
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> _______________________________________________ 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>