Skip to contents

Hurdle Negative Binomial Family

Source: R/family_zi.R
ratiod_hurdle_negbin.Rd

Hurdle model where zero vs positive counts are modelled as separate processes.

$$P(Y = 0) = 1 - \theta$$ $$P(Y = y | Y > 0) = P_{NB}^+(y), \quad y > 0$$

where \(\theta\) is the probability of a positive count and \(P_{NB}^+\) is the truncated-at-zero negative binomial.

Usage

ratiod_hurdle_negbin(
  link_num = "log",
  link_denom = "log",
  link_hurdle = "logit",
  denom_family = "negbin"
)

Arguments

Link function for the count mean: "log".

Link function for the denominator mean: "log" (default); a gamma denominator also takes "identity" or "inverse".

Link function for the hurdle 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.

Value

A ratiod_family object with hurdle structure

Details

The hurdle model differs from zero-inflation in interpretation:

Zero-inflated: Zeros come from two sources (structural + sampling)

  • \(\pi\): probability of structural zero

  • Some zeros also come from the count process

Hurdle: Zeros are all from one source (the "hurdle" process)

  • \(1 - \theta\): probability of zero

  • All positive counts from truncated count process

Choose hurdle when zeros represent a distinct biological/sampling state.

Examples

# Create family object
fam <- ratiod_hurdle_negbin()
print(fam)
#> tulpaRatio family: hurdle_negbin_negbin 
#> [Hurdle model]
#> Hurdle negative binomial numerator, negative binomial denominator 
#> 
#> Numerator:  neg_binomial_2(log)
#> Hurdle:     bernoulli(logit)
#> Denominator: neg_binomial_2 (log) 

# Simulate hurdle data (presence/absence + abundance)
set.seed(123)
n <- 60
presence_prob <- 0.7
df <- data.frame(
  count = ifelse(runif(n) > presence_prob, 0,
                 rnbinom(n, size = 3, mu = 10) + 1),
  total = rnbinom(n, size = 5, mu = 50),
  habitat = factor(rep(c("forest", "grassland"), each = n/2)),
  site = factor(rep(1:10, each = n/10))
)

# \donttest{
# Fit model
fit <- tratio(
  count | total ~ habitat + (1 | site),
  data = df,
  family = ratiod_hurdle_negbin(),
  control = list(iter = 400, chains = 1, seed = 1)
)
summary(fit)
#> ratio model summary
#> ===================
#> 
#> Inference: Exact (Tier 1) via hmc
#> Family: hurdle_negbin_negbin 
#> 
#> Fixed effects (numerator):
#>             parameter   mean    sd   2.5% 97.5%  rhat ess_bulk
#>       num:(Intercept)  2.273 0.145  2.003 2.528 0.998       65
#>  num:habitatgrassland -0.015 0.196 -0.418 0.348 1.007       60
#> 
#> Fixed effects (denominator):
#>               parameter   mean    sd   2.5% 97.5%  rhat ess_bulk
#>       denom:(Intercept)  3.939 0.126  3.723 4.190 1.003      106
#>  denom:habitatgrassland -0.001 0.157 -0.310 0.287 1.001      106
#> 
#> Variance components:
#>  parameter  mean    sd  2.5% 97.5%  rhat ess_bulk
#>  phi_denom 3.854 0.635 2.724 5.169 1.009      144
#>    phi_num 4.744 1.466 2.562 7.899 1.020      144
#>   sigma_re 0.092 0.071 0.008 0.292 1.009       50
#> 
#> Zero/one inflation:
#>   parameter  mean    sd 2.5% 97.5%  rhat ess_bulk
#>  beta_zi[1] 0.904 0.249 0.43 1.428 1.001      123
#> 
#> Diagnostics:
#>   Divergences: 0 
#>   Avg. acceptance: 0.851 
#>   Chains: 1 
#>   Iterations: 400 
#>   Warmup: 200 
#>   R-hat > 1.01 for 1 parameter(s), largest 1.020 (phi_num).
#>   Bulk ESS < 100 for 3 parameter(s), smallest 50 (sigma_re).
#>   Run longer chains or reparameterise before using these estimates.
# }