Zero-Inflation and Hurdle Models

Installation

install.packages("pak")                   # development version
pak::pak("gcol33/tulpaRatio")

Overview

Ecological and survey data often exhibit excess zeros beyond what standard distributions predict. tulpaRatio provides several families to handle this:

Zero-inflated models assume zeros come from two sources: 1. Structural zeros (species truly absent, observer non-response) 2. Sampling zeros (from the count/binomial process)

Hurdle models assume all zeros come from a single “hurdle” process, with positive counts from a truncated distribution.

Available Families

Count Data (Numerator)

Family	Base	ZI Type	Use Case
`ratiod_zinegbin()`	Negative binomial	Mixture	Overdispersed counts with excess zeros
`ratiod_zipois()`	Poisson	Mixture	Counts where overdispersion is from zeros only
`ratiod_hurdle_negbin()`	Negative binomial	Hurdle	Presence/absence + abundance
`ratiod_hurdle_pois()`	Poisson	Hurdle	Presence/absence + Poisson counts

Proportion Data (Binomial Numerator)

Family	Base	ZI Type	Use Case
`ratiod_zibinomial()`	Binomial	Zero-inflated	Excess zero success rates
`ratiod_oibinomial()`	Binomial	One-inflated	Excess 100% success rates
`ratiod_zoibinomial()`	Binomial	Zero-and-one	Excess at both boundaries
`ratiod_hurdle_binomial()`	Binomial	Hurdle	Zero vs positive successes

Zero-Inflated Binomial

Use ratiod_zibinomial() when modeling proportions with excess zeros beyond what binomial predicts.

Example: Species Detection

library(tulpaRatio)

ctl <- list(seed = 1)

# Simulate detection data with structural zeros
set.seed(123)
n <- 60
habitat <- factor(rep(c("forest", "grassland"), each = n/2))
site <- factor(rep(1:10, each = n/10))
site_eff <- rnorm(10, 0, 0.5)[site]
# Detection probability rises in grassland and varies by site
p <- plogis(-0.4 + 0.5 * (habitat == "grassland") + site_eff)
# Structural zeros: 15% in forest, 40% in grassland
zi_prob <- ifelse(habitat == "forest", 0.15, 0.4)

detection <- data.frame(
  detected = ifelse(
    runif(n) < zi_prob,
    0,  # Structural zero
    rbinom(n, size = 10, prob = p)  # Binomial detections
  ),
  trials = 10,
  habitat = habitat,
  site = site
)

# Fit zero-inflated binomial
fit_zi <- tratio(
  detected | trials ~ habitat + (1 | site),
  zi = zi_binomial(~ habitat),  # Predictors for zero-inflation
  data = detection,
  family = ratiod_zibinomial(),
  control = ctl
)

print(coef(summary(fit_zi)), digits = 3)

The zi argument takes a specification such as zi_binomial(), whose formula is the linear predictor for the zero-inflation probability (beta_zi above, on the logit scale). Use zi_binomial(~ 1) for intercept-only or include covariates that predict structural zeros.

One-Inflated Binomial

Use ratiod_oibinomial() when data has excess ones (100% success rates).

Example: Survey Completion

# Survey completion data with some "always complete" respondents
set.seed(456)
n <- 60
oi_prob <- 0.2  # 20% always complete all items

survey <- data.frame(
  completed = ifelse(
    runif(n) < oi_prob,
    20,  # Structural one (all items completed)
    rbinom(n, size = 20, prob = 0.7)
  ),
  total_items = 20,
  age_group = factor(sample(c("young", "middle", "old"), n, replace = TRUE)),
  education = rnorm(n)
)

# Fit one-inflated binomial
fit_oi <- tratio(
  completed | total_items ~ age_group + education,
  data = survey,
  family = ratiod_oibinomial(),
  control = ctl
)

print(coef(summary(fit_oi)), digits = 3)

Zero-and-One Inflated Binomial (ZOIB)

Use ratiod_zoibinomial() when data has excess at both boundaries.

Example: Vegetation Cover

Vegetation cover data often has structural zeros (no vegetation) and structural ones (complete cover).

set.seed(789)
n <- 60
zi_prob <- 0.15  # 15% bare ground
oi_prob <- 0.10  # 10% complete cover (given not bare)

cover <- data.frame(
  cover_points = sapply(1:n, function(i) {
    if (runif(1) < zi_prob) {
      0  # Structural zero
    } else if (runif(1) < oi_prob) {
      100  # Structural one
    } else {
      rbinom(1, size = 100, prob = 0.5)
    }
  }),
  total_points = 100,
  elevation = rnorm(n),  # standardized elevation
  aspect = factor(sample(c("N", "S", "E", "W"), n, replace = TRUE))
)

# Fit ZOIB model
fit_zoib <- tratio(
  cover_points | total_points ~ elevation + aspect,
  zi = zoib(
    zi_formula = ~ elevation,  # Zero-inflation predictors
    oi_formula = ~ elevation   # One-inflation predictors
  ),
  data = cover,
  family = ratiod_zoibinomial(),
  control = ctl
)

print(coef(summary(fit_zoib)), digits = 3)

Hurdle Binomial

Use ratiod_hurdle_binomial() when zeros represent a distinct process from positive counts.

Example: Occupancy with Detection

set.seed(321)
n <- 60
occupancy <- data.frame(
  visits = 5L,
  habitat_quality = rnorm(n),
  site = factor(rep(1:12, each = 5))
)
# Presence rises with habitat quality; where present, the detections out of
# five visits are a binomial count given at least one detection, with a
# site effect on the detection probability
present <- runif(n) < plogis(0.8 + 0.8 * occupancy$habitat_quality)
p_detect <- plogis(0.3 + 0.4 * occupancy$habitat_quality +
                     rnorm(12, 0, 0.5)[occupancy$site])
detect_positive <- function(p) {
  repeat {
    y <- rbinom(1, 5, p)
    if (y > 0) return(y)
  }
}
occupancy$detections <- ifelse(present, vapply(p_detect, detect_positive, numeric(1)), 0)

# Fit hurdle binomial
fit_hurdle <- tratio(
  detections | visits ~ habitat_quality + (1 | site),
  zi = hurdle_binomial(~ habitat_quality),  # Predictors for hurdle (presence)
  data = occupancy,
  family = ratiod_hurdle_binomial(),
  control = ctl
)

print(coef(summary(fit_hurdle)), digits = 3)

Zero-Inflated Count Models

For count numerators with continuous denominators (CPUE-type data):

# Zero-inflated CPUE data
set.seed(111)
n <- 60
season <- factor(rep(c("spring", "summer"), each = n/2))
vessel <- factor(rep(1:12, each = n/12))
# A vessel effect shared by catch and effort
vessel_eff <- rnorm(12, 0, 0.5)[vessel]
depth <- rnorm(n)
# Structural zeros: 15% in spring, 35% in summer
zi_prob <- ifelse(season == "spring", 0.15, 0.35)

cpue <- data.frame(
  catch = ifelse(
    runif(n) < zi_prob,
    0,
    rpois(n, exp(2 - 0.3 * depth + vessel_eff))
  ),
  effort = rgamma(n, shape = 4, rate = exp(-vessel_eff)),
  depth = depth,
  season = season,
  vessel = vessel
)

# Zero-inflated Poisson with gamma denominator
fit_zi_cpue <- tratio(
  catch | effort ~ depth + season + (1 | vessel),
  zi = zi_poisson(~ season),
  data = cpue,
  family = ratiod_zipois(denom_family = "gamma"),
  control = ctl
)

print(coef(summary(fit_zi_cpue)), digits = 3)

Choosing Between Models

Scenario	Recommended Family
Zeros from mixed sources (structural + sampling)	`ratiod_zi*()`
All zeros from single process (presence/absence)	`ratiod_hurdle_*()`
Excess at 0 only	`ratiod_zibinomial()`
Excess at 1 (100%) only	`ratiod_oibinomial()`
Excess at both 0 and 1	`ratiod_zoibinomial()`
Overdispersion primarily from zeros	`ratiod_zipois()`
Overdispersion in zeros AND positive counts	`ratiod_zinegbin()`

Model Comparison

Compare ZI and non-ZI models using LOO-CV:

# Fit standard model
fit_standard <- tratio(
  detected | trials ~ habitat + (1 | site),
  data = detection,
  family = ratiod_binomial(),
  control = ctl
)

# Fit zero-inflated model
fit_zi <- tratio(
  detected | trials ~ habitat + (1 | site),
  zi = zi_binomial(~ 1),
  data = detection,
  family = ratiod_zibinomial(),
  control = ctl
)

# Compare
ratiod_compare(fit_standard, fit_zi)

The Pareto \(k\) warning marks observations whose site effect moves too far when they are left out for importance sampling to follow; their contributions to the comparison are less reliable than the rest.

Posterior Predictive Checks

Validate ZI models with posterior predictive checks:

# Check proportion of zeros
pp_check(fit_zi, type = "stat", stat = function(y) mean(y == 0))

# Check distribution
pp_check(fit_zi, type = "dens_overlay")