R/zero-inflation.R

In DynCount: Bayesian Dynamic Models for Count Time Series

# Zero-inflation diagnostics --------------------------------------------------

#' Posterior probability that each observed zero is structural
#'
#' For a zero-inflated fit (Poisson or binomial), returns for every observed
#' zero the posterior probability that it is a *structural* zero (gate closed)
#' rather than an ordinary *sampling* zero generated by the observation process.
#' A single gate-open probability governs all observations (it does not vary
#' over time or with covariates), while the gate itself is drawn separately for
#' every observation.
#'
#' @details
#' The model introduces a latent gate indicator \eqn{v_t \in \{0, 1\}}: when
#' the gate is open (\eqn{v_t = 1}) the count comes from the observation model
#' (Poisson or binomial); when closed (\eqn{v_t = 0}) the count is a structural
#' zero. The sampler stores \eqn{v_t} for each draw, so the structural-zero
#' probability is
#' \deqn{\Pr(\text{structural} \mid y_t = 0) = 1 - \overline{v_t},}
#' the posterior mean of the gate being closed. Non-zero observations are
#' always sampling observations and have structural probability `0`.
#'
#' @param object A `"dynamic_fit"` object fitted with `zeros = "inflated"`
#'   (Poisson or binomial family).
#' @param zeros_only If `TRUE` (default), return only the rows where the
#'   observation is zero; with `FALSE`, return one row per observation (the
#'   non-zero ones with `p_structural = 0`).
#'
#' @return A data frame with columns `time`, `observed`, `p_structural`
#'   (posterior probability the zero is structural) and `p_sampling`
#'   (`1 - p_structural`).
#'
#' @examples
#' sim <- simulate_dynamic_poisson(60, 0.2, 2, zero_inflation = 0.25, seed = 1)
#' fit <- fit_dynamic_model(sim$y, zero_inflation = TRUE, nsave = 300, nburn = 200,
#'                          seed = 1)
#' structural_zero_prob(fit)
#' @export
structural_zero_prob <- function(object, zeros_only = TRUE) {
  if (!inherits(object, "dynamic_fit")) stop("`object` must be a dynamic_fit.", call. = FALSE)
  if (identical(object$spec$family, "multinomial")) {
    stop("structural_zero_prob() is not available for the multinomial family ",
         "(zero inflation is not supported).", call. = FALSE)
  }
  if (object$spec$zeros != "inflated") {
    stop("structural_zero_prob() requires a model fitted with `zeros = \"inflated\"`.",
         call. = FALSE)
  }
  gate_open <- colMeans(object$draws$gate)        # P(gate open | data)
  y <- object$data$y
  p_struct <- ifelse(y == 0, 1 - gate_open, 0)
  out <- data.frame(
    time = seq_along(y),
    observed = y,
    p_structural = p_struct,
    p_sampling = 1 - p_struct
  )
  if (zeros_only) out <- out[y == 0, , drop = FALSE]
  rownames(out) <- NULL
  out
}