[R-sig-ME] Calculating Confidence intervals for glmmTMB predictions

[Michael Whitby]

I am trying to calculate confidence intervals of predictions made from a
glmmTMB model with zero inflation using an ar1 covariance structure. The
model uses complex bases (both poly and scale) in both the model and zi
formula.

I have looked through a few issues posted on github and the original paper
describing glmmTMB. However, many methods are proposed and shown, and it
has only caused me confustion. Additionally, it appears glmmTMB has changed
slightly (e.g. ziform depricated, and addition of allow.new.levels
argurment), leading me to think that the correct method may be different.

Here is what I've looked at so far:
https://github.com/glmmTMB/glmmTMB/issues/324
https://github.com/glmmTMB/glmmTMB/issues/378
https://www.biorxiv.org/content/biorxiv/suppl/2017/05/01/132753.DC1/132753-2.pdf
https://r-sig-mixed-models.r-project.narkive.com/EOC1ab9c/r-sig-me-zero-inflated-glmmtmb-with-poly-confidence-band

Currently, I am simply taking the prediction and se on the response scale
(since this should include the ZI term as well), converting back to the log
scale, calculating the 0.95 confidence interval and converting back to the
response scale.

Here is an example using the Salamanders data set (this is somewhat
simplified as it doesn't use ar1 or complex basis, but I think it still
demonstrates my question). I would like some reassurance that I calculated
the upper and lower limits correctly.

library(glmmTMB)
data("Salamanders")

# View(Salamanders)

fit = glmmTMB(count~spp + cover + mined + poly(DOP, 3)+(1|site),
              ziformula=~spp + mined,
              dispformula = ~DOY,
              data = Salamanders,
              family=nbinom2)
# DOP
-----------------------------------------------------------------------------------------


newdata2 <- expand.grid(
  site = "new",
  spp = unique(Salamanders$spp),
  mined = factor(c("yes", "no"), levels = levels(Salamanders$mined)),
  cover = mean(Salamanders$cover),
  DOY = mean(Salamanders$DOY),
  DOP = seq(from = min(Salamanders$DOP), to=max(Salamanders$DOP), length =
25)
)

preds2 <- predict(fit, newdata2, se=T, allow.new.levels = T,
type='response')
newdata2$pred=preds2$fit
newdata2$se = preds2$se.fit
newdata2$ulimit = exp(log(newdata2$pred)+(qnorm(0.975)*log(newdata2$se)))
newdata2$llimit = exp(log(newdata2$pred)-(qnorm(0.975)*log(newdata2$se)))

library(ggplot2)

ggplot(data=newdata2, aes(x=DOP, y = pred) )+
  geom_ribbon(aes(ymin=llimit, ymax=ulimit, fill=mined), alpha = 0.25)+
  geom_line(aes(color = mined), size=1)+
  facet_wrap(~spp)



Michael Whitby
michael.whitby using gmail.com
609-923-0973

	[[alternative HTML version deleted]]