Calculating Confidence intervals for glmmTMB predictions

Hi Michael,

as the zero-inflation component and count component of the model work in
"opposite directions" (a higher expected value for the zero inflation means
a lower response, but a higher value for the conditional model means a
higher response), the recommendation from the package authors (see
especially
https://github.com/glmmTMB/glmmTMB/issues/378#issuecomment-417850022, which
refers to the Appendix A in the paper) is to use a bootstrapping-method to
compute the CI (see
https://journal.r-project.org/archive/2017/RJ-2017-066/RJ-2017-066.pdf, page
391-392).

I have implemented various options lately in my ggeffects-package
(https://strengejacke.github.io/ggeffects). There are examples for the
various predict-options (count component only, conditioning on
zero-inflation, taking random effects into account etc.) in the package
vignettes, see here:
https://strengejacke.github.io/ggeffects/articles/randomeffects.html#margina
l-effects-for-zero-inflated-mixed-models

These features are currently only available in the GitHub-version of
ggeffects, but a CRAN-submission is planned in the next days. ggpredict()
will do all the work for you and allows you to easily create a plot. A
function-call would look like this:

ggpredict(fit, c("DOP", "mined", "spp"), type = "fe.zi") %>% plot()


If you want to this this manually, you could use this code:

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)
)

prdat.sim <- get_glmmTMB_predictions(fit, newdata2, 1000) sims <-
exp(prdat.sim$cond) * (1 - stats::plogis(prdat.sim$zi))

newdata2$predicted <- apply(sims, 1, mean) newdata2$conf.low <- apply(sims,
1, quantile, probs = .025) newdata2$conf.high <- apply(sims, 1, quantile,
probs = .975) newdata2$std.error <- apply(sims, 1, sd)

ggplot(data=newdata2, aes(x=DOP, y = predicted) )+
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=mined), alpha = 0.25)+
  geom_line(aes(color = mined), size=1)+
  facet_wrap(~spp)

get_glmmTMB_predictions <- function(model, newdata, nsim) {
  tryCatch(
    {
      x.cond <- stats::model.matrix(lme4::nobars(stats::formula(model)[-2]),
newdata)
      beta.cond <- lme4::fixef(model)$cond
      
      ziformula <- model$modelInfo$allForm$ziformula
      x.zi <- stats::model.matrix(lme4::nobars(stats::formula(ziformula)),
newdata)
      beta.zi <- lme4::fixef(model)$zi
      
      pred.condpar.psim <- MASS::mvrnorm(n = nsim, mu = beta.cond, Sigma =
stats::vcov(model)$cond)
      pred.cond.psim <- x.cond %*% t(pred.condpar.psim)
      pred.zipar.psim <- MASS::mvrnorm(n = nsim, mu = beta.zi, Sigma =
stats::vcov(model)$zi)
      pred.zi.psim <- x.zi %*% t(pred.zipar.psim)
      
      list(cond = pred.cond.psim, zi = pred.zi.psim)
    },
    error = function(x) { NULL },
    warning = function(x) { NULL },
    finally = function(x) { NULL }
  )
}


Best
Daniel

-----Urspr?ngliche Nachricht-----
Von: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> Im
Auftrag von Michael Whitby
Gesendet: Dienstag, 1. Januar 2019 19:37
An: r-sig-mixed-models at r-project.org
Betreff: [R-sig-ME] Calculating Confidence intervals for glmmTMB predictions

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-infl
ated-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 at gmail.com
609-923-0973


_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

--

_____________________________________________________________________

Universit?tsklinikum Hamburg-Eppendorf; K?rperschaft des ?ffentlichen Rechts; Gerichtsstand: Hamburg | www.uke.de
Vorstandsmitglieder: Prof. Dr. Burkhard G?ke (Vorsitzender), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Pr?l?, Marya Verdel
_____________________________________________________________________

SAVE PAPER - THINK BEFORE PRINTING