I am trying to implement a zero inflated model in INLA.
I know a basic zero inflated Poisson can be implemented with "zeroinflatedpoisson1" as the family argument as described here (where I also took the example below). However, this paper describes a different approach where the two likelihoods are written into the INLA call, and the response is an matrix with $2n$ rows and 2 columns corresponding to the likelihoods.
I have tried to implement it below with the Owls data from glmmTMB, but the results suggest it has not implemented two likelihoods. However, the results are close, but different, from the Poisson version shown at the end.
library(INLA)
library(tidyverse)
df <- glmmTMB::Owls |> mutate(NCalls= SiblingNegotiation)
df1 <- df |> mutate(ziy = ifelse(NCalls == 0, 0 , 1),
NCalls = NA)
df$ziy <- NA
dff <- rbind(df1,df)
dff$Y <- data.frame(dff[, c("ziy", "NCalls")])
# zip model - doesnt seem to work
m.zip <- inla(Y~ FoodTreatment + ArrivalTime + SexParent +
logBroodSize + f(Nest, model = "iid"),
family= c("binomial", "poisson"),
data=dff,
control.fixed=list(mean=0, prec = 1),
control.predictor=list(compute=TRUE),
control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE))
summary(m.zip)
# Call:
# c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data,
# quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials
# = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates,
# verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute =
# control.compute, ", " control.predictor = control.predictor, control.family = control.family,
# ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode =
# control.mode, control.expert = control.expert, ", " control.hazard = control.hazard,
# control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale =
# control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
# ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", "
# blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", "
# silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
# .parent.frame)")
# Time used:
# Pre = 2.82, Running = 0.42, Post = 0.0163, Total = 3.26
# Fixed effects:
# mean sd 0.025quant 0.5quant 0.975quant mode kld
# (Intercept) 4.549 0.404 3.759 4.547 5.350 4.544 0
# FoodTreatmentSatiated -0.634 0.036 -0.705 -0.634 -0.564 -0.634 0
# ArrivalTime -0.134 0.009 -0.153 -0.134 -0.116 -0.134 0
# SexParentMale 0.049 0.036 -0.022 0.049 0.119 0.049 0
# logBroodSize 0.518 0.238 0.041 0.521 0.981 0.526 0
#
# Random effects:
# Name Model
# Nest IID model
#
# Model hyperparameters:
# mean sd 0.025quant 0.5quant 0.975quant mode
# Precision for Nest 4.75 1.60 2.24 4.54 8.46 4.13
#
# Deviance Information Criterion (DIC) ...............: 5564.07
# Deviance Information Criterion (DIC, saturated) ....: 3883.28
# Effective number of parameters .....................: 28.46
#
# Watanabe-Akaike information criterion (WAIC) ...: 5685.81
# Effective number of parameters .................: 140.35
#
# Marginal log-Likelihood: -2833.00
# CPO, PIT is computed
# Posterior summaries for the linear predictor and the fitted values are computed
# (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
# Poisson
m.p <- inla(NCalls ~ FoodTreatment + ArrivalTime + SexParent +
logBroodSize + f(Nest, model = "iid"),
family= c("poisson"),
data=df,
control.fixed=list(mean=0, prec = 1),
control.predictor=list(compute=TRUE),
control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE))
summary(m.p)
# Call:
# c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data,
# quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials
# = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates,
# verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute =
# control.compute, ", " control.predictor = control.predictor, control.family = control.family,
# ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode =
# control.mode, control.expert = control.expert, ", " control.hazard = control.hazard,
# control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale =
# control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
# ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", "
# blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", "
# silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
# .parent.frame)")
# Time used:
# Pre = 2.86, Running = 0.379, Post = 0.0433, Total = 3.28
# Fixed effects:
# mean sd 0.025quant 0.5quant 0.975quant mode kld
# (Intercept) 4.457 0.383 3.707 4.455 5.215 4.452 0
# FoodTreatmentSatiated -0.585 0.036 -0.656 -0.585 -0.515 -0.585 0
# ArrivalTime -0.130 0.009 -0.148 -0.130 -0.111 -0.130 0
# SexParentMale 0.042 0.036 -0.029 0.042 0.113 0.042 0
# logBroodSize 0.511 0.219 0.072 0.513 0.936 0.518 0
#
# Random effects:
# Name Model
# Nest IID model
#
# Model hyperparameters:
# mean sd 0.025quant 0.5quant 0.975quant mode
# Precision for Nest 5.92 2.01 2.77 5.65 10.57 5.14
#
# Deviance Information Criterion (DIC) ...............: 4945.64
# Deviance Information Criterion (DIC, saturated) ....: 3264.82
# Effective number of parameters .....................: 28.40
#
# Watanabe-Akaike information criterion (WAIC) ...: 5060.26
# Effective number of parameters .................: 135.53
#
# Marginal log-Likelihood: -2520.86
# CPO, PIT is computed
# Posterior summaries for the linear predictor and the fitted values are computed
# (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
#
