example of simulation ("parametric Bayes")-based confidence intervals for Z-I gamma

zigamma_pred.R

library(glmmTMB)
set.seed(101)
dd <- data.frame(x = rnorm(1000))

dd$y0 <- simulate_new( ~ 1,
                     seed = 101,
                     ziformula = ~ 1,
                     family = ziGamma(link = "log"),
                     newdata = dd,
                     newparams = list(
                         beta = c(0),
                         betadisp = 1,
                         betazi = 0
                     )
                     )[[1]]

all.equal(plogis(0), mean(dd$y0 == 0), tolerance = 0.02)
par(las=1); hist(dd$y0[dd$y0 > 0], breaks = 50, freq = FALSE)
curve(dgamma(x, shape = exp(1), scale = 1/exp(1)), add = TRUE, col = "red", lwd = 2)

## is zigamma simulation working right?
params <- list(
    beta = c(0,1),
    betazi = c(-2,1),
    betadisp = 1)

dd$y <- simulate_new( ~ 1 + x,
                     seed = 101,
                     ziformula = ~ 1 + x,
                     family = ziGamma(link = "log"),
                     newdata = dd,
                     newparams = params
                     )[[1]]
dd <- dd[order(dd$y), ]
plot(dd$y)

m <- glmmTMB(y ~ x, ziformula = ~ x,
             family = ziGamma(link = "log"),
             data = dd)

mm <- unlist(fixef(m))
cbind(mm, true = unlist(params))

vv <- vcov(m, full = TRUE)

pars <- MASS::mvrnorm(1000, mm, vv)
## surprisingly slow (2 minutes)
system.time(
    simres <- apply(pars, 1,
                    \(p) predict(m, newparams = p, type = "response")
                    )
)

X <- getME(m, "X")
Xzi <- getME(m, "Xzi")
mypred <- function(p) {
    zprob <- plogis(Xzi %*% p[3:4])
    mu <- exp(X %*% p[1:2])
    mu * (1-zprob)
}
simres2 <- apply(pars, 1, mypred)

plot(mypred(pars[100,]))

lines(dd$y, col = 2, lwd = 2)

## matplot(simres, pch = ".", col = "black")
matplot(simres2, pch = ".", col = "black")
matpoints(simres, pch = ".", col = "red")
points(dd$y, col = 2, pch = 16)
points(predict(m, type = "response"), col = adjustcolor("blue", alpha = 0.2),
       pch = 16)