Zero-inflated negative binomial for the numerator process. Models excess zeros as a mixture:
$$P(Y = 0) = \pi + (1 - \pi) \cdot P_{NB}(0)$$ $$P(Y = y) = (1 - \pi) \cdot P_{NB}(y), \quad y > 0$$
where \(\pi\) is the zero-inflation probability and \(P_{NB}\) is the negative binomial PMF.
Usage
ratiod_zinegbin(
link_num = "log",
link_denom = "log",
link_zi = "logit",
denom_family = "negbin"
)Arguments
- link_num
Link function for the count mean: "log".
- link_denom
Link function for the denominator mean: "log" (default); a gamma denominator also takes "identity" or "inverse".
- link_zi
Link function for the zero-inflation probability: "logit" (default), "probit" or "cloglog".
- denom_family
Denominator distribution: "negbin" (default), "gamma" for continuous effort, or "fixed" for a known exposure. A known exposure is data: it scales the count mean, which is then the rate times the exposure, and the ratio is the rate.
Details
The zero-inflation probability can have its own linear predictor via the
zi argument in tratio(). If not specified, a single intercept is estimated.
For ecological data, zero-inflation often represents:
Structural zeros (species truly absent vs not detected)
Sampling zeros (inadequate effort)
False negatives
See also
ratiod_hurdle_negbin() for hurdle model alternative
Examples
# Create family object
fam <- ratiod_zinegbin()
print(fam)
#> tulpaRatio family: zinegbin_negbin
#> [Zero-inflated model]
#> Zero-inflated negative binomial numerator, negative binomial denominator
#>
#> Numerator: neg_binomial_2(log)
#> ZI prob: bernoulli(logit)
#> Denominator: neg_binomial_2 (log)
# Simulate zero-inflated count data
set.seed(123)
n <- 60
zi_prob <- 0.3
habitat <- factor(rep(c("forest", "grassland"), each = n/2))
site <- factor(rep(1:10, each = n/10))
u <- rnorm(10, sd = 0.5)[site]
df <- data.frame(
count = ifelse(runif(n) < zi_prob, 0,
rnbinom(n, size = 3, mu = exp(2 + 0.5 * (habitat == "grassland") + u))),
total = rnbinom(n, size = 5, mu = 50 * exp(u)),
habitat = habitat,
site = site
)
# \donttest{
# Fit model
fit <- tratio(
count | total ~ habitat + (1 | site),
data = df,
family = ratiod_zinegbin(),
control = list(iter = 400, chains = 1, seed = 1)
)
summary(fit)
#> ratio model summary
#> ===================
#>
#> Inference: Exact (Tier 1) via hmc
#> Family: zinegbin_negbin
#>
#> Fixed effects (numerator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> num:(Intercept) 1.918 0.298 1.471 2.587 1.013 65
#> num:habitatgrassland 0.341 0.419 -0.523 1.206 1.004 59
#>
#> Fixed effects (denominator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> denom:(Intercept) 4.080 0.263 3.679 4.586 1.049 51
#> denom:habitatgrassland -0.289 0.399 -1.164 0.561 1.007 48
#>
#> Variance components:
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> phi_denom 3.496 0.655 2.316 4.963 1.008 91
#> phi_num 3.126 1.118 0.945 5.543 1.000 52
#> sigma_re 0.622 0.173 0.342 0.956 1.050 37
#>
#> Zero/one inflation:
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> beta_zi[1] -0.855 0.337 -1.553 -0.315 1.025 39
#>
#> Diagnostics:
#> Divergences: 0
#> Avg. acceptance: 0.928
#> Chains: 1
#> Iterations: 400
#> Warmup: 200
#> R-hat > 1.01 for 4 parameter(s), largest 1.050 (sigma_re).
#> Bulk ESS < 100 for 8 parameter(s), smallest 37 (sigma_re).
#> Run longer chains or reparameterise before using these estimates.
# }