Add zero-truncated gamma-Poisson functions

[sebdalgarno]

Summary

• Adds log_lik_gamma_pois_zt(), dev_gamma_pois_zt(), ran_gamma_pois_zt(), and res_gamma_pois_zt() for the zero-truncated gamma-Poisson (negative binomial) distribution, mirroring the existing _zi family
• Likelihood is the standard gamma-Poisson density divided by 1 - P(0; lambda, theta)
• Random samples use rejection sampling
• Residuals centre on the truncated mean lambda / (1 - P(0)); standardised form scales by the analytically-derived truncated variance
• Deviances reuse the untruncated NB deviance approximation since the saturated MLE for the truncated likelihood has no closed form

Why

Useful for fitting zero-truncated negative binomial models (e.g., count data filtered to be >= 1 by design) where current extras lacks residual / random-draw / log-likelihood support. Encountered in a kelp biomass analysis where plant frond counts are filtered to fronds >= 1.

Test plan

• Unit tests added in test-log-lik.R, test-dev.R, test-ran.R, test-res.R following _zi patterns
• All targeted tests pass (devtools::test(filter = "log-lik|dev|ran|res"))
• Documentation regenerated (devtools::document())
• NEWS.md updated
• Verified rejection sampling guarantees Y >= 1
• Verified standardised residuals on simulated data approximate N(0, 1)

🤖 Generated with Claude Code

[Copilot]

Pull request overview

Adds support for the zero-truncated gamma-Poisson (negative binomial) distribution to extras, mirroring the existing gamma-Poisson family helpers.

Changes:

• Added new _zt functions: log_lik_gamma_pois_zt(), dev_gamma_pois_zt(), ran_gamma_pois_zt(), res_gamma_pois_zt().
• Added unit tests covering log-likelihood, deviance, random draws, and residual behavior for the new _zt family.
• Updated exports and regenerated documentation to include the new functions (plus a NEWS entry).

Reviewed changes

Copilot reviewed 57 out of 57 changed files in this pull request and generated 8 comments.

Show a summary per file

File	Description
R/dev.R	Adds dev_gamma_pois_zt() implementation.
R/log-lik.R	Adds log_lik_gamma_pois_zt() implementation.
R/ran.R	Adds ran_gamma_pois_zt() via rejection sampling.
R/res.R	Adds res_gamma_pois_zt() with multiple residual types and simulation option.
NAMESPACE	Exports the new _zt functions.
NEWS.md	Adds a user-facing NEWS entry for the new functions.
tests/testthat/test-dev.R	Adds tests for _zt deviance and residual sign behavior.
tests/testthat/test-log-lik.R	Adds tests for _zt log-likelihood correctness and vectorization.
tests/testthat/test-ran.R	Adds tests for _zt sampling behavior and truncated mean.
tests/testthat/test-res.R	Adds tests for _zt residual types and simulated standardized residuals.
man/dev_gamma_pois_zi.Rd	Cross-links deviance docs to include dev_gamma_pois_zt().
man/dev_gamma_pois_zt.Rd	New generated documentation for dev_gamma_pois_zt().
man/log_lik_bern.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_beta.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_beta_binom.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_binom.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_exp.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_gamma.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_gamma_pois.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_gamma_pois_zi.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_gamma_pois_zt.Rd	New generated documentation for log_lik_gamma_pois_zt().
man/log_lik_lnorm.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_neg_binom.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_norm.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_pois.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_pois_zi.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_skewnorm.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_student.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/log_lik_unif.Rd	Cross-links log-likelihood docs to include log_lik_gamma_pois_zt().
man/ran_bern.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_beta_binom.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_binom.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_gamma.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_gamma_pois.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_gamma_pois_zi.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_gamma_pois_zt.Rd	New generated documentation for ran_gamma_pois_zt().
man/ran_lnorm.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_neg_binom.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_norm.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_pois.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_pois_zi.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_skewnorm.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/ran_student.Rd	Cross-links random-sampling docs to include ran_gamma_pois_zt().
man/res_bern.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_beta_binom.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_binom.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_gamma.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_gamma_pois.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_gamma_pois_zi.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_gamma_pois_zt.Rd	New generated documentation for res_gamma_pois_zt().
man/res_lnorm.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_neg_binom.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_norm.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_pois.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_pois_zi.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_skewnorm.Rd	Cross-links residual docs to include res_gamma_pois_zt().
man/res_student.Rd	Cross-links residual docs to include res_gamma_pois_zt().

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[sebdalgarno]

prrof from claude on why the MLE approach is better than previous approach:

#
# Closed-form `dev_nb(x, lambda, theta) + 2 * (log(1 - p0(x)) - log(1 - p0(lambda)))`
# treats lambda_sat = x as the saturated point, but under zero-truncation
# the saturated MLE is not x. The result is deviance residuals with
# variance != 1 even when the model is correctly specified.
#
# Numeric saturated MLE produces deviance residuals with variance ~ 1 as
# expected, and matches the analytic standardized residual variance.

library(extras)

# Closed-form (incorrect under truncation)
dev_zt_old <- function(x, lambda, theta) {
  dev_nb <- extras::dev_neg_binom(x, lambda = lambda, theta = theta, res = FALSE)
  log_p0_mu <- (-1 / theta) * log1p(lambda * theta)
  log_p0_y  <- (-1 / theta) * log1p(x * theta)
  dev_nb + 2 * (log1p(-exp(log_p0_y)) - log1p(-exp(log_p0_mu)))
}

# Numeric saturated MLE (correct)
log_lik_zt <- function(x, lambda, theta) {
  ll <- stats::dnbinom(x, mu = lambda, size = 1 / theta, log = TRUE)
  l0 <- stats::dnbinom(0, mu = lambda, size = 1 / theta, log = TRUE)
  ll - log1p(-exp(l0))
}
dev_zt_new <- function(x, lambda, theta) {
  upper <- max(x)
  opt_lambda <- extras:::parallel_optimize(
    f = function(lam) -log_lik_zt(x, lambda = lam, theta = theta),
    interval = c(.Machine$double.eps, upper),
    N = length(x)
  )
  d <- log_lik_zt(x, opt_lambda, theta) - log_lik_zt(x, lambda, theta)
  d[!is.na(d) & d < 0 & d > -1e-7] <- 0
  2 * d
}

# Sample ZT-NB
ran_zt <- function(n, lambda, theta) {
  out <- integer(n); rem <- seq_len(n)
  while (length(rem)) {
    d <- as.integer(rnbinom(length(rem), mu = lambda, size = 1 / theta))
    a <- d > 0; out[rem[a]] <- d[a]; rem <- rem[!a]
  }
  out
}

# Analytic standardized-residual variance under truncation -- should be 1
trunc_stats <- function(lambda, theta) {
  log_p0 <- (-1 / theta) * log1p(lambda * theta)
  one_minus_p0 <- -expm1(log_p0)
  mu_t <- lambda / one_minus_p0
  var_t <- (lambda + theta * lambda^2 + lambda^2) / one_minus_p0 - mu_t^2
  list(mu = mu_t, var = var_t)
}

# 30 samples of size n at true (lambda, theta); model is correctly specified
set.seed(1)
nrep <- 30; n <- 4000; lambda <- 2.5; theta <- 0.6
ts <- trunc_stats(lambda, theta)

v_old <- v_new <- v_std <- numeric(nrep)
for (k in seq_len(nrep)) {
  x <- ran_zt(n, lambda, theta)
  v_old[k] <- var(sign(x - ts$mu) * sqrt(pmax(dev_zt_old(x, lambda, theta), 0)))
  v_new[k] <- var(sign(x - ts$mu) * sqrt(pmax(dev_zt_new(x, lambda, theta), 0)))
  v_std[k] <- var((x - ts$mu) / sqrt(ts$var))
}

cat(sprintf("variance of deviance residuals over %d ZT-NB samples (n = %d, model correctly specified):\n", nrep, n))
cat(sprintf("  closed-form:    mean = %.3f  sd = %.3f\n", mean(v_old), sd(v_old)))
cat(sprintf("  numeric MLE:    mean = %.3f  sd = %.3f\n", mean(v_new), sd(v_new)))
cat(sprintf("  standardized:   mean = %.3f  sd = %.3f  (analytic check; should be ~1)\n",
            mean(v_std), sd(v_std)))