---
title: "Zero-Inflation and Hurdle Models"
author: "Gilles Colling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Zero-Inflation and Hurdle Models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  eval = identical(Sys.getenv("NOT_CRAN"), "true") || identical(Sys.getenv("IN_PKGDOWN"), "true")
)
```

## Installation

```{r, eval = FALSE}
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

```{r}
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

```{r}
# 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).

```{r}
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

```{r}
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):

```{r}
# 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:

```{r}
# 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:

```{r}
# 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")
```


## See Also

- `vignette("getting-started")` - Basic usage
- `vignette("spatial-temporal")` - Spatial and temporal models
- `?ratiod_zibinomial`, `?ratiod_zoibinomial` - Function documentation