The V Lab
AudienceLearners in public health, epidemiology, medicine, and health data science
Study timeApproximately 150–210 minutes
PrerequisitesBasic regression, likelihood, confidence intervals, and R

About the simulated data This module uses a fixed-seed simulated public-health cohort of respiratory emergency-department visits. Follow-up differs in person-years, and some people have structural zeros because their care is outside the participating-capture system. The data contain no real people or clinical records.

How to use this tutorial

Define the event and the exposure → inspect the count distribution → fit a Poisson rate model → assess dispersion and zeros → choose an interpretable extension → check predictions and calibration → report the design and limitations.

The word “zero-inflated” names a probability model, not a diagnosis. Start with the data-generating and measurement processes; then ask whether a mixture model answers a meaningful question.

Learning objectives

After completing this tutorial, you should be able to:

  • distinguish a count, a rate, a risk, and a probability of zero;
  • specify and interpret a Poisson regression with offset(log(person_years));
  • convert a log-rate coefficient to an incidence rate ratio (IRR), confidence interval, and absolute prediction;
  • calculate Pearson dispersion and explain common sources of overdispersion;
  • compare Poisson, quasi-Poisson, negative binomial, and zero-inflated Poisson (ZIP) approaches;
  • write the ZIP likelihood and interpret both its count and structural-zero components;
  • obtain ZIP predictions for the susceptible count mean, structural-zero probability, marginal mean, and observed-zero probability;
  • assess observed-versus-predicted frequency distributions and calibration without treating AIC as a scientific verdict;
  • distinguish ZIP from a hurdle model;
  • recognize the limits imposed by missing data, clustering, study design, and causal questions; and
  • produce an auditable count-model report.

1 The question, outcome, and denominator

1.1 Count, rate, risk, and probability are different quantities

A count is the number of observed events: 3 respiratory emergency-department (ED) visits. A rate divides events by person-time: 3 visits in 1.5 person-years equals 2 visits per person-year. A risk is the probability of at least one event in a specified interval, which needs a fixed follow-up definition and is not generally the same as a rate. A zero probability, P(Y=0)P(Y=0), is neither an IRR nor a risk unless the outcome and observation window make it one.

For person ii, let YiY_i be ED visits and tit_i be observed person-years. A Poisson rate model assumes

Yi∣Xi∼Poisson⁡(μi),log⁡(μi)=log⁡(ti)+β0+XiTβ. Y_i\mid X_i \sim \operatorname{Poisson}(\mu_i), \qquad \log(\mu_i)=\log(t_i)+\beta_0+X_i^T\beta.

The coefficient of log(t_i) is fixed at one. It is an offset, not a covariate whose association is estimated. Thus exp⁡(XiTβ)\exp(X_i^T\beta) is an expected event rate per person-year, and exponentiated covariate coefficients are IRRs holding other modeled predictors fixed.

Do not divide the outcome and then also use an offset For unequal observation time, model the integer count with offset(log(person_years)). A linear model for visits/person_years has the wrong error structure, gives unstable weight to short follow-up, and should not be combined with the offset as though it were a second correction.

1.2 The teaching cohort

The simulated cohort has n=2400n=2400 people. Program participation, smoking, severity, rural residence, and age influence the rate while access barriers and rural residence influence whether care is outside capture. This mechanism deliberately creates more zeros than a one-process Poisson model expects.

summary_table <- data.frame(
  Quantity = c("People", "Total visits", "Person-years", "Observed zero visits",
               "Mean visits", "Variance of visits", "Mean follow-up"),
  Value = c(nrow(count_data), sum(count_data$visits), sum(count_data$person_years),
            sum(count_data$visits == 0), mean(count_data$visits), var(count_data$visits),
            mean(count_data$person_years))
)
knitr::kable(summary_table, digits = 2,
  caption = "Descriptive summary of the simulated respiratory ED-visit cohort")
Descriptive summary of the simulated respiratory ED-visit cohort
Quantity Value
People 2400.00
Total visits 1717.00
Person-years 2926.23
Observed zero visits 1451.00
Mean visits 0.72
Variance of visits 1.29
Mean follow-up 1.22 
count_tab <- table(factor(count_data$visits, levels = 0:max(count_data$visits)))
barplot(count_tab, col = palette_count["teal"], border = "white", xlab = "Observed ED visits",
        ylab = "Number of people", main = "Observed count distribution")
A histogram of nonnegative ED-visit counts shows that zero is the most common value, followed by progressively fewer observations at one, two, and higher counts.

The distribution of observed respiratory ED-visit counts has a large bar at zero and a decreasing right tail. This plot motivates diagnosis but does not by itself identify a zero-inflated process.

1.3 Before modeling: data audit and complete cases

The outcome must be a nonnegative integer; identify whether a zero means no event, no opportunity, a recording rule, or an unobserved event. Confirm that person-time is positive and belongs to the same observation period as the count. Check duplicate rows, impossible dates, coding changes, and whether care outside a network is observed.

analysis_variables <- c("visits", "person_years", "program", "smoking", "severity",
                        "rural", "age_c10", "access_barrier")
complete_rows <- complete.cases(count_data[, analysis_variables])
data.frame(
  Eligible_rows = nrow(count_data), Complete_rows = sum(complete_rows),
  Excluded_for_missingness = sum(!complete_rows),
  Minimum_person_years = min(count_data$person_years),
  Maximum_person_years = max(count_data$person_years)
)
Eligible_rows Complete_rows Excluded_for_missingness Minimum_person_years Maximum_person_years
2400 2400 0 0.45 2

This simulated example is complete by construction. In a study, complete-case analysis is unbiased only under restrictive conditions; report missingness by variable and use a prespecified approach such as multiple imputation when justified. Do not impute the event count or person-time casually: the imputation model must honor their measurement process.

2 Poisson regression for rates

2.1 Fit the model and inspect the linear predictor

poisson_fit <- glm(
  visits ~ program + smoking + severity + rural + age_c10 + offset(log(person_years)),
  family = poisson(link = "log"), data = count_data
)
poisson_fit
## 
## Call:  glm(formula = visits ~ program + smoking + severity + rural + 
##     age_c10 + offset(log(person_years)), family = poisson(link = "log"), 
##     data = count_data)
## 
## Coefficients:
##    (Intercept)  programProgram      smokingYes        severity      ruralRural  
##         -0.588          -0.254           0.371           0.368          -0.176  
##        age_c10  
##          0.124  
## 
## Degrees of Freedom: 2399 Total (i.e. Null);  2394 Residual
## Null Deviance:       3490 
## Residual Deviance: 3110  AIC: 5400

The program coefficient compares the adjusted visit rate in the program group with the standard group. The simulated mechanism includes a program log-rate effect of −0.30-0.30, but real analyses never have this privileged access to the data-generating truth.

2.2 IRRs, confidence intervals, and absolute predictions

For coefficient β̂j\hat\beta_j, the IRR is exp⁡(β̂j)\exp(\hat\beta_j). A Wald 95% confidence interval is exp⁡{β̂j±1.96SE⁡(β̂j)}\exp\{\hat\beta_j \pm 1.96\operatorname{SE}(\hat\beta_j)\}. This is a conditional relative comparison, not an individual probability and not necessarily a causal effect.

poisson_ci <- confint.default(poisson_fit)
poisson_irr <- data.frame(
  Term = names(coef(poisson_fit)), IRR = exp(coef(poisson_fit)),
  `Lower 95% CI` = exp(poisson_ci[, 1]), `Upper 95% CI` = exp(poisson_ci[, 2]),
  check.names = FALSE
)
poisson_effects <- poisson_irr[poisson_irr$Term != "(Intercept)", ]
knitr::kable(poisson_effects, digits = 3,
  caption = "Poisson-model covariate incidence rate ratios with Wald 95% confidence intervals")
Poisson-model covariate incidence rate ratios with Wald 95% confidence intervals
Term IRR Lower 95% CI Upper 95% CI
programProgram programProgram 0.776 0.705 0.853
smokingYes smokingYes 1.449 1.308 1.605
severity severity 1.444 1.377 1.515
ruralRural ruralRural 0.839 0.757 0.929
age_c10 age_c10 1.132 1.087 1.180

The intercept is omitted from the IRR table. Its exponentiated value is the baseline rate for the reference levels when continuous covariates equal zero, not a covariate-effect ratio.

reference_profiles <- data.frame(
  person_years = 1,
  program = factor(c("Standard", "Program"), levels = levels(count_data$program)),
  smoking = factor("No", levels = levels(count_data$smoking)), severity = 0,
  rural = factor("Urban", levels = levels(count_data$rural)), age_c10 = 0,
  access_barrier = factor("No", levels = levels(count_data$access_barrier))
)
reference_profiles$poisson_mean_visits <- predict(poisson_fit, reference_profiles, type = "response")
knitr::kable(reference_profiles[, c("program", "person_years", "poisson_mean_visits")], digits = 3,
  caption = "Poisson-model expected visits for two reference profiles with one person-year of follow-up")
Poisson-model expected visits for two reference profiles with one person-year of follow-up
program person_years poisson_mean_visits
Standard 1 0.555
Program 1 0.431

Absolute predictions make the baseline, follow-up time, and covariate pattern visible. For a Poisson model, the predicted probability of at least one visit is 1−exp⁡(−μ̂)1-\exp(-\hat\mu); it is not the same quantity as the expected count μ̂\hat\mu.

3 Dispersion and model extensions

3.1 Pearson dispersion is a diagnostic, not a ritual

Under a correctly specified Poisson model, conditional variance equals conditional mean. A common descriptive statistic is

ϕ̂=∑irPi2n−p,rPi=yi−μ̂iμ̂i. \hat\phi=\frac{\sum_i r_{Pi}^2}{n-p}, \qquad r_{Pi}=\frac{y_i-\hat\mu_i}{\sqrt{\hat\mu_i}}.

Values well above one suggest overdispersion, but a single cutoff cannot determine its cause. Omitted predictors, nonlinear effects, heterogeneous rates, repeated observations, dependence within clinics, recording mixtures, and excess zeros can all inflate residual variation.

pearson_dispersion <- sum(residuals(poisson_fit, type = "pearson")^2) / df.residual(poisson_fit)
data.frame(
  Pearson_dispersion = pearson_dispersion,
  Residual_degrees_of_freedom = df.residual(poisson_fit),
  Interpretation = ifelse(pearson_dispersion > 1.2, "Overdispersion is present; investigate its source.",
                          "No large dispersion signal in this descriptive check.")
)
Pearson_dispersion Residual_degrees_of_freedom Interpretation
1.36 2394 Overdispersion is present; investigate its source.

3.2 Quasi-Poisson and negative binomial models

A quasi-Poisson model retains the mean model E(Y∣X)=μE(Y\mid X)=\mu but estimates a scale factor for standard errors. It has no full likelihood, so ordinary AIC comparisons are unavailable. A negative binomial (NB) model adds a dispersion parameter, often written Var⁡(Y∣X)=μ+μ2/θ\operatorname{Var}(Y\mid X)=\mu+\mu^2/\theta. It is especially useful when unobserved multiplicative heterogeneity is a plausible rate mechanism.

quasi_fit <- glm(
  visits ~ program + smoking + severity + rural + age_c10 + offset(log(person_years)),
  family = quasipoisson(link = "log"), data = count_data
)
if (!requireNamespace("MASS", quietly = TRUE)) {
  stop("This module requires MASS for the negative binomial example. Install MASS first.", call. = FALSE)
}
nb_fit <- MASS::glm.nb(
  visits ~ program + smoking + severity + rural + age_c10 + offset(log(person_years)),
  data = count_data
)
comparison_rates <- data.frame(
  Model = c("Poisson", "Quasi-Poisson", "Negative binomial"),
  Program_IRR = c(exp(coef(poisson_fit)["programProgram"]),
                  exp(coef(quasi_fit)["programProgram"]), exp(coef(nb_fit)["programProgram"])),
  Program_SE_log_IRR = c(summary(poisson_fit)$coefficients["programProgram", "Std. Error"],
                        summary(quasi_fit)$coefficients["programProgram", "Std. Error"],
                        summary(nb_fit)$coefficients["programProgram", "Std. Error"])
)
knitr::kable(comparison_rates, digits = 3,
  caption = "Program rate-ratio estimates from models that address variance differently")
Program rate-ratio estimates from models that address variance differently
Model Program_IRR Program_SE_log_IRR
Poisson 0.776 0.049
Quasi-Poisson 0.776 0.057
Negative binomial 0.779 0.061

An NB model is not a zero-inflated model NB allows greater count variance through continuous rate heterogeneity. ZIP adds a separate point mass for structural zeros. Both may fit the same data, neither is automatically the truth, and their coefficients answer different model-based questions.

4 Zero-inflated Poisson models

4.1 When a two-process story is plausible

In the simulated setting, an observed zero can arise because a person is outside participating-care capture (a structural zero) or because a person inside capture has no visit during follow-up (a Poisson zero). That is a substantive mixture story. Excess observed zeros alone do not prove that story: a wrong mean model, short follow-up, unmodeled clustering, or NB heterogeneity can produce the same visual symptom.

Let Si=1S_i=1 identify the structural-zero state and let πi=P(Si=1∣Zi)\pi_i=P(S_i=1\mid Z_i). Conditional on Si=0S_i=0, Yi∼Poisson⁡(μi)Y_i\sim\operatorname{Poisson}(\mu_i). The ZIP probability mass function is

P(Yi=0)=πi+(1−πi)e−μi,P(Yi=y>0)=(1−πi)e−μiμiyy!. P(Y_i=0)=\pi_i+(1-\pi_i)e^{-\mu_i}, \qquad P(Y_i=y>0)=(1-\pi_i)\frac{e^{-\mu_i}\mu_i^y}{y!}.

We use log links for the count mean and a logit link for πi\pi_i. The count component describes the susceptible/captured process; the zero component describes membership in the always-zero state. A positive zero-component coefficient increases structural-zero odds, not the odds of any observed zero.

4.2 Fit a transparent ZIP likelihood

The following functions implement maximum likelihood directly with optim(). For zeros, the likelihood is evaluated as log⁡[π+(1−π)e−μ]\log[\pi+(1-\pi)e^{-\mu}], avoiding unstable multiplication on the ordinary probability scale. The inverse Hessian provides approximate standard errors, so convergence and the validity of the Hessian must be checked.

# This teaching ZIP implementation uses direct maximum likelihood, stable
# log-probabilities, multiple starts, and Hessian-based standard errors.
fit_zip <- function(count_formula, zero_formula, data, maxit = 5000) {
  count_frame <- model.frame(count_formula, data = data, na.action = na.fail)
  count_terms <- terms(count_frame)
  y <- model.response(count_frame)
  if (any(y < 0 | y != floor(y))) stop("The outcome must contain nonnegative integer counts.")
  X <- model.matrix(count_terms, count_frame)
  count_offset <- model.offset(count_frame)
  if (is.null(count_offset)) count_offset <- rep(0, length(y))
  zero_frame <- model.frame(zero_formula, data = data, na.action = na.fail)
  zero_terms <- terms(zero_frame)
  Z <- model.matrix(zero_terms, zero_frame)
  if (nrow(Z) != length(y)) stop("Count and zero formulas produced different rows.")
  poisson_start <- glm(count_formula, family = poisson(link = "log"), data = data)
  beta_start <- coef(poisson_start)
  poisson_zero <- mean(exp(-fitted(poisson_start)))
  pi_start <- (mean(y == 0) - poisson_zero) / (1 - poisson_zero)
  pi_start <- pmin(pmax(pi_start, 0.02), 0.80)
  n_count <- ncol(X)
  logspace_add <- function(a, b) {
    larger <- pmax(a, b)
    larger + log(exp(a - larger) + exp(b - larger))
  }
  nll <- function(par) {
    beta <- par[seq_len(n_count)]
    gamma <- par[n_count + seq_len(ncol(Z))]
    eta_count <- drop(count_offset + X %*% beta)
    if (any(!is.finite(eta_count)) || any(abs(eta_count) > 700)) {
      return(.Machine$double.xmax / 100)
    }
    mu <- exp(eta_count)
    eta_zero <- drop(Z %*% gamma)
    log_pi <- plogis(eta_zero, log.p = TRUE)
    log_one_minus_pi <- plogis(-eta_zero, log.p = TRUE)
    is_zero <- y == 0
    log_lik <- numeric(length(y))
    log_lik[is_zero] <- logspace_add(
      log_pi[is_zero], log_one_minus_pi[is_zero] - mu[is_zero]
    )
    log_lik[!is_zero] <- log_one_minus_pi[!is_zero] +
      dpois(y[!is_zero], lambda = mu[!is_zero], log = TRUE)
    if (any(!is.finite(log_lik))) return(.Machine$double.xmax / 100)
    -sum(log_lik)
  }

  zero_intercept_starts <- unique(c(qlogis(pi_start), qlogis(c(0.05, 0.25, 0.50))))
  starts <- lapply(
    zero_intercept_starts,
    function(intercept) c(beta_start, intercept, rep(0, ncol(Z) - 1L))
  )
  candidates <- lapply(starts, function(start) {
    optim(
      par = start, fn = nll, method = "BFGS", hessian = TRUE,
      control = list(maxit = maxit, reltol = 1e-10)
    )
  })
  valid <- vapply(
    candidates,
    function(candidate) candidate$convergence == 0 && is.finite(candidate$value),
    logical(1)
  )
  if (!any(valid)) stop("No ZIP optimization start converged.", call. = FALSE)
  valid_candidates <- candidates[valid]
  opt <- valid_candidates[[which.min(vapply(valid_candidates, `[[`, numeric(1), "value"))]]

  hessian <- (opt$hessian + t(opt$hessian)) / 2
  hessian_eigen <- eigen(hessian, symmetric = TRUE, only.values = TRUE)$values
  hessian_condition <- max(hessian_eigen) / min(hessian_eigen)
  if (any(!is.finite(hessian_eigen)) || min(hessian_eigen) <= 0 ||
      !is.finite(hessian_condition) || hessian_condition > 1e10) {
    stop("The ZIP Hessian is not sufficiently positive definite.", call. = FALSE)
  }
  vcov_matrix <- solve(hessian)

  numerical_gradient <- vapply(seq_along(opt$par), function(j) {
    step <- 1e-5 * max(1, abs(opt$par[j]))
    upper <- lower <- opt$par
    upper[j] <- upper[j] + step
    lower[j] <- lower[j] - step
    (nll(upper) - nll(lower)) / (2 * step)
  }, numeric(1))
  gradient_max <- max(abs(numerical_gradient))
  if (!is.finite(gradient_max) || gradient_max > 0.01) {
    stop("The ZIP solution has an unacceptably large numerical gradient.", call. = FALSE)
  }
  coefficient_names <- c(paste0("count_", colnames(X)), paste0("zero_", colnames(Z)))
  names(opt$par) <- coefficient_names
  dimnames(vcov_matrix) <- list(coefficient_names, coefficient_names)
  factor_levels <- function(frame) lapply(frame[vapply(frame, is.factor, logical(1))], levels)
  structure(list(coefficients = opt$par, vcov = vcov_matrix, logLik = -opt$value,
    convergence = opt$convergence, count_terms = count_terms, zero_terms = zero_terms,
    count_xlevels = factor_levels(count_frame), zero_xlevels = factor_levels(zero_frame),
    count_contrasts = attr(X, "contrasts"), zero_contrasts = attr(Z, "contrasts"),
    n_count = n_count, n_zero = ncol(Z), nobs = length(y), y = y,
    starts_tried = length(starts), gradient_max = gradient_max,
    hessian_min_eigen = min(hessian_eigen), hessian_condition = hessian_condition),
    class = "teaching_zip")
}

predict_zip <- function(object, newdata) {
  count_terms <- delete.response(object$count_terms)
  count_frame <- model.frame(count_terms, data = newdata, na.action = na.pass,
                             xlev = object$count_xlevels)
  X <- model.matrix(count_terms, count_frame, contrasts.arg = object$count_contrasts)
  count_offset <- model.offset(count_frame)
  if (is.null(count_offset)) count_offset <- rep(0, nrow(X))
  zero_frame <- model.frame(object$zero_terms, data = newdata, na.action = na.pass,
                            xlev = object$zero_xlevels)
  Z <- model.matrix(object$zero_terms, zero_frame, contrasts.arg = object$zero_contrasts)
  beta <- object$coefficients[seq_len(object$n_count)]
  gamma <- object$coefficients[object$n_count + seq_len(object$n_zero)]
  mu <- exp(drop(count_offset + X %*% beta))
  pi <- plogis(drop(Z %*% gamma))
  data.frame(count_mean = mu, structural_zero_probability = pi,
             marginal_mean = (1 - pi) * mu,
             observed_zero_probability = pi + (1 - pi) * exp(-mu))
}
zip_fit <- fit_zip(
  visits ~ program + smoking + severity + rural + age_c10 + offset(log(person_years)),
  ~ access_barrier + rural, data = count_data
)
if (zip_fit$convergence != 0) stop("ZIP optimization did not converge.")
zip_diagnostics <- data.frame(
  Quantity = c("Convergence code", "Starting values tried", "Maximum absolute gradient",
               "Minimum Hessian eigenvalue", "Hessian condition number"),
  Value = c(zip_fit$convergence, zip_fit$starts_tried, zip_fit$gradient_max,
            zip_fit$hessian_min_eigen, zip_fit$hessian_condition)
)
knitr::kable(zip_diagnostics, digits = 4,
  caption = "ZIP optimization and curvature diagnostics")
ZIP optimization and curvature diagnostics
Quantity Value
Convergence code 0.0000
Starting values tried 4.0000
Maximum absolute gradient 0.0004
Minimum Hessian eigenvalue 17.7536
Hessian condition number 170.9350 
zip_se <- sqrt(diag(zip_fit$vcov))
zip_ci <- cbind(zip_fit$coefficients - 1.96 * zip_se, zip_fit$coefficients + 1.96 * zip_se)
zip_table <- data.frame(
  Component = ifelse(grepl("^count_", names(zip_fit$coefficients)), "Count mean", "Structural-zero probability"),
  Term = sub("^(count_|zero_)", "", names(zip_fit$coefficients)),
  Estimate = zip_fit$coefficients, SE = zip_se,
  `Exp(estimate)` = exp(zip_fit$coefficients),
  `Lower 95% CI` = exp(zip_ci[, 1]), `Upper 95% CI` = exp(zip_ci[, 2]), check.names = FALSE
)
zip_effect_table <- zip_table[zip_table$Term != "(Intercept)", ]
knitr::kable(zip_effect_table, digits = 3,
  caption = "ZIP covariate effects. Exponentiated count coefficients are IRRs; exponentiated zero-component coefficients are structural-zero odds ratios.")
ZIP covariate effects. Exponentiated count coefficients are IRRs; exponentiated zero-component coefficients are structural-zero odds ratios.
Component Term Estimate SE Exp(estimate) Lower 95% CI Upper 95% CI
count_programProgram Count mean programProgram -0.262 0.053 0.77 0.693 0.854
count_smokingYes Count mean smokingYes 0.355 0.058 1.43 1.273 1.596
count_severity Count mean severity 0.375 0.027 1.46 1.380 1.534
count_ruralRural Count mean ruralRural 0.148 0.067 1.16 1.017 1.322
count_age_c10 Count mean age_c10 0.096 0.023 1.10 1.053 1.152
zero_access_barrierYes Structural-zero probability access_barrierYes 1.330 0.159 3.78 2.768 5.167
zero_ruralRural Structural-zero probability ruralRural 0.752 0.181 2.12 1.488 3.023

The omitted count intercept exponentiates to the reference profile’s baseline rate per person-year; the omitted zero-component intercept exponentiates to the reference profile’s baseline structural-zero odds. Neither intercept is a covariate-effect ratio.

The programProgram count coefficient estimates a program IRR within the modeled nonstructural process, conditional on covariates. The access-barrier coefficient in the zero component changes odds of being structurally zero. Neither automatically estimates the population-average causal impact of program participation.

4.3 Predictions need all four quantities

For a covariate profile, ZIP returns:

  1. μ\mu: mean count among the nonstructural process;
  2. π\pi: probability of structural zero;
  3. (1−π)μ(1-\pi)\mu: marginal expected observed count; and
  4. π+(1−π)e−μ\pi+(1-\pi)e^{-\mu}: predicted observed-zero probability.
zip_profiles <- reference_profiles
zip_profiles <- zip_profiles[rep(seq_len(nrow(zip_profiles)), each = 2), ]
rownames(zip_profiles) <- NULL
zip_profiles$access_barrier <- factor(
  rep(c("No", "Yes"), times = 2),
  levels = levels(count_data$access_barrier)
)
zip_prediction <- cbind(zip_profiles[, c("program", "access_barrier", "person_years")],
                        predict_zip(zip_fit, zip_profiles))
knitr::kable(zip_prediction, digits = 3,
  caption = "ZIP predictions for reference profiles. The marginal mean and observed-zero probability combine both model components.")
ZIP predictions for reference profiles. The marginal mean and observed-zero probability combine both model components.
program access_barrier person_years count_mean structural_zero_probability marginal_mean observed_zero_probability
Standard No 1 0.750 0.187 0.610 0.571
Standard Yes 1 0.750 0.465 0.402 0.717
Program No 1 0.577 0.187 0.470 0.643
Program Yes 1 0.577 0.465 0.309 0.765 
# Parametric coefficient simulation preserves covariance between both ZIP components.
set.seed(20260821)
n_draws <- 1200
coefficient_draws <- matrix(
  rnorm(n_draws * length(zip_fit$coefficients)), nrow = n_draws
) %*% chol(zip_fit$vcov)
coefficient_draws <- sweep(coefficient_draws, 2, zip_fit$coefficients, "+")

count_terms_profiles <- delete.response(zip_fit$count_terms)
count_frame_profiles <- model.frame(
  count_terms_profiles, zip_profiles, xlev = zip_fit$count_xlevels
)
X_profiles <- model.matrix(
  count_terms_profiles, count_frame_profiles,
  contrasts.arg = zip_fit$count_contrasts
)
profile_offset <- model.offset(count_frame_profiles)
Z_profiles <- model.matrix(
  zip_fit$zero_terms,
  model.frame(zip_fit$zero_terms, zip_profiles, xlev = zip_fit$zero_xlevels),
  contrasts.arg = zip_fit$zero_contrasts
)

beta_draws <- coefficient_draws[, seq_len(zip_fit$n_count), drop = FALSE]
gamma_draws <- coefficient_draws[, zip_fit$n_count + seq_len(zip_fit$n_zero), drop = FALSE]
mu_draws <- exp(sweep(beta_draws %*% t(X_profiles), 2, profile_offset, "+"))
pi_draws <- plogis(gamma_draws %*% t(Z_profiles))
marginal_draws <- (1 - pi_draws) * mu_draws
zero_probability_draws <- pi_draws + (1 - pi_draws) * exp(-mu_draws)

zip_prediction[["Marginal mean lower 95% limit"]] <- apply(marginal_draws, 2, quantile, 0.025)
zip_prediction[["Marginal mean upper 95% limit"]] <- apply(marginal_draws, 2, quantile, 0.975)
zip_prediction[["Zero probability lower 95% limit"]] <- apply(zero_probability_draws, 2, quantile, 0.025)
zip_prediction[["Zero probability upper 95% limit"]] <- apply(zero_probability_draws, 2, quantile, 0.975)
knitr::kable(zip_prediction, digits = 3,
  caption = "ZIP profile predictions with approximate 95% coefficient-uncertainty intervals")
ZIP profile predictions with approximate 95% coefficient-uncertainty intervals
program access_barrier person_years count_mean structural_zero_probability marginal_mean observed_zero_probability Marginal mean lower 95% limit Marginal mean upper 95% limit Zero probability lower 95% limit Zero probability upper 95% limit
Standard No 1 0.750 0.187 0.610 0.571 0.555 0.668 0.543 0.600
Standard Yes 1 0.750 0.465 0.402 0.717 0.347 0.462 0.680 0.755
Program No 1 0.577 0.187 0.470 0.643 0.427 0.515 0.616 0.668
Program Yes 1 0.577 0.465 0.309 0.765 0.267 0.357 0.733 0.796

These intervals propagate model-coefficient uncertainty while retaining the cross-covariance between the count and zero components. They are not prediction intervals for a future person’s realized count and do not include model-selection or structural uncertainty.

Do not report μ\mu as the expected count for everyone In ZIP, μ\mu is the Poisson mean within the nonstructural component. The population-average expected observed count is (1−π)μ(1-\pi)\mu. Clearly state which quantity your prediction table reports.

5 Diagnostics, calibration, and comparisons

5.1 Frequency distribution: observed versus predicted

We compare the observed count frequencies with model-based expected frequencies. For a Poisson model, P(Y=k)=e−μμk/k!P(Y=k)=e^{-\mu}\mu^k/k!. For ZIP, use the mixture formula for zero and multiply ordinary Poisson probabilities by (1−π)(1-\pi) for positive counts. The last category pools the tail to make a readable diagnostic.

zip_pred_all <- predict_zip(zip_fit, count_data)
mu_pois <- fitted(poisson_fit)
max_show <- 5
observed_frequency <- tabulate(
  pmin(count_data$visits, max_show + 1) + 1,
  nbins = max_show + 2
)
poisson_frequency <- sapply(0:max_show, function(k) sum(dpois(k, mu_pois)))
poisson_frequency <- c(poisson_frequency, nrow(count_data) - sum(poisson_frequency))
zip_frequency <- sapply(0:max_show, function(k) {
  if (k == 0) sum(zip_pred_all$observed_zero_probability) else {
    sum((1 - zip_pred_all$structural_zero_probability) * dpois(k, zip_pred_all$count_mean))
  }
})
zip_frequency <- c(zip_frequency, nrow(count_data) - sum(zip_frequency))
frequency_check <- data.frame(
  Count_category = c(as.character(0:max_show), paste0(max_show + 1, "+")),
  Observed = observed_frequency, Poisson_expected = poisson_frequency, ZIP_expected = zip_frequency
)
knitr::kable(frequency_check, digits = 1,
  caption = "Observed and model-predicted count frequencies, with the upper tail pooled")
Observed and model-predicted count frequencies, with the upper tail pooled
Count_category Observed Poisson_expected ZIP_expected
0 1451 1276.1 1445.9
1 498 718.8 506.2
2 257 274.0 257.8
3 120 90.9 113.2
4 47 28.3 46.6
5 15 8.5 18.6
6+ 12 3.4 11.7 
barplot(
  t(as.matrix(frequency_check[, c("Observed", "Poisson_expected", "ZIP_expected")])),
  beside = TRUE, names.arg = frequency_check$Count_category,
  col = c(palette_count["navy"], palette_count["orange"], palette_count["teal"]),
  border = NA, xlab = "Visit-count category", ylab = "Frequency",
  main = "Observed versus expected frequencies"
)
legend("topright", legend = c("Observed", "Poisson expected", "ZIP expected"),
       fill = c(palette_count["navy"], palette_count["orange"], palette_count["teal"]), bty = "n")
A grouped bar chart shows observed, Poisson-expected, and ZIP-expected frequencies for zero through five or more ED visits. The Poisson model underpredicts zero visits, while the ZIP model is closer across the displayed categories.

Observed count frequencies are compared with frequencies expected under the Poisson and ZIP models. Agreement across zero, small positive counts, and the pooled tail is more informative than a single fit statistic.

5.2 Calibration is a prediction question

Calibration asks whether predicted outcomes agree with observed outcomes on the relevant scale and population. Group people by predicted marginal mean or zero probability, then compare observed and expected values. Internal apparent calibration can look excellent even when a model is overfit; use resampling or an external cohort for prediction claims.

calibration_group <- cut(rank(zip_pred_all$marginal_mean, ties.method = "first"),
                         breaks = quantile(seq_len(nrow(count_data)), probs = seq(0, 1, 0.1)),
                         include.lowest = TRUE, labels = FALSE)
calibration <- aggregate(cbind(observed_visits = count_data$visits,
                               predicted_visits = zip_pred_all$marginal_mean,
                               observed_zero = count_data$visits == 0,
                               predicted_zero = zip_pred_all$observed_zero_probability) ~ calibration_group,
                         FUN = mean)
knitr::kable(calibration, digits = 3,
  caption = "Apparent ZIP calibration by deciles of predicted marginal mean")
Apparent ZIP calibration by deciles of predicted marginal mean
calibration_group observed_visits predicted_visits observed_zero predicted_zero
1 0.179 0.179 0.850 0.851
2 0.300 0.281 0.767 0.779
3 0.392 0.361 0.717 0.731
4 0.396 0.450 0.717 0.682
5 0.546 0.536 0.646 0.641
6 0.600 0.639 0.613 0.587
7 0.729 0.768 0.546 0.534
8 0.987 0.943 0.458 0.470
9 1.167 1.189 0.417 0.421
10 1.858 1.811 0.317 0.329 
plot(calibration$predicted_visits, calibration$observed_visits, pch = 19,
     col = palette_count["teal"], xlab = "Mean predicted visits", ylab = "Mean observed visits",
     main = "Apparent calibration of the ZIP marginal mean")
abline(0, 1, lty = 2, col = palette_count["gray"], lwd = 2)
A scatterplot compares average predicted ZIP marginal mean on the horizontal axis with average observed ED visits on the vertical axis for ten prediction groups. A diagonal reference line indicates perfect calibration, and the points lie near it.

Calibration of the ZIP marginal mean across ten groups ranked by prediction. Points near the diagonal indicate agreement between average predicted and observed visit counts within these simulated data.

5.3 AIC helps compare likelihoods, not theories

Poisson, NB, and ZIP likelihood-based AIC values can be compared only when they are fitted to the same outcome, rows, and likelihood convention. Quasi-Poisson has no ordinary likelihood AIC. A lower AIC reflects an estimated out-of-sample information criterion, not proof that structural zeros exist, that causal confounding is solved, or that predictions transport to another health system.

zip_aic <- -2 * zip_fit$logLik + 2 * length(zip_fit$coefficients)
aic_table <- data.frame(
  Model = c("Poisson", "Negative binomial", "ZIP"),
  AIC = c(AIC(poisson_fit), AIC(nb_fit), zip_aic),
  Note = c("Single Poisson count process", "Extra-Poisson rate heterogeneity", "Poisson process plus structural-zero mixture")
)
knitr::kable(aic_table, digits = 1,
  caption = "AIC values for likelihood-based models fitted to the same simulated cohort")
AIC values for likelihood-based models fitted to the same simulated cohort
Model AIC Note
Poisson 5400 Single Poisson count process
Negative binomial 5234 Extra-Poisson rate heterogeneity
ZIP 5086 Poisson process plus structural-zero mixture

6 ZIP versus hurdle models

A hurdle model also has two components, but all zeros come from the binary hurdle and the positive-count distribution is zero-truncated. Once the hurdle is crossed, a positive count is guaranteed. ZIP instead allows a susceptible person to have zero events through the Poisson process. The distinction should follow the process: a screening rule that blocks all recorded events may support a hurdle; participating care with a chance of no visit supports a ZIP-like mixture.

Question ZIP Hurdle model
Can the count process generate zero? Yes No; it is truncated at zero
What does the binary part model? Structural-zero membership Zero versus any positive count
What does a positive count imply? Not structural zero The hurdle was crossed
Is excess zeros enough to choose it? No No

7 Design, clustering, and causal boundaries

7.1 The regression coefficient is usually an association

The program variable is randomized in this simulation only for illustration. In an observational study, program participation may reflect severity, access, prior utilization, clinic policy, or variables not measured. Adjustment can describe conditional associations; a causal rate effect needs an explicit target trial, temporal ordering, a well-defined intervention, no important unmeasured confounding, appropriate handling of selection and censoring, positivity, and a valid outcome-capture process.

Repeated visits within people, people within clinics, or neighborhoods within regions violate independence. Cluster-robust standard errors can address some inference issues with enough independent clusters; mixed-effects or generalized estimating equation approaches address different estimands and correlation structures. A ZIP likelihood that assumes independent rows is not repaired merely by adding many covariates.

Measurement systems can create the zeros Claims records, network boundaries, insurance changes, distance, and access barriers can determine whether an event is visible. Treat zero patterns as possible information-system behavior, not as intrinsic patient behavior. This matters for equity, transportability, and causal interpretation.

8 A complete analysis and reporting template

8.1 Summarize the main simulated analysis

program_irr_pois <- exp(coef(poisson_fit)["programProgram"])
program_irr_nb <- exp(coef(nb_fit)["programProgram"])
program_irr_zip <- exp(zip_fit$coefficients["count_programProgram"])
case_summary <- data.frame(
  Result = c("Participants", "Observed zero-visit proportion", "Pearson dispersion", "Poisson program IRR", "NB program IRR", "ZIP count-component program IRR", "ZIP convergence"),
  Value = c(nrow(count_data), mean(count_data$visits == 0), pearson_dispersion,
            program_irr_pois, program_irr_nb, program_irr_zip, zip_fit$convergence)
)
knitr::kable(case_summary, digits = 3,
  caption = "Selected results from the simulated respiratory ED-visit analysis")
Selected results from the simulated respiratory ED-visit analysis
Result Value
Participants 2400.000
Observed zero-visit proportion 0.605
Pearson dispersion 1.363
Poisson program IRR 0.776
NB program IRR 0.779
ZIP count-component program IRR 0.770
ZIP convergence 0.000

In this simulated cohort, 60.5% of people had no recorded respiratory ED visit. The Poisson Pearson dispersion was 1.36, indicating that a one-process Poisson variance was inadequate. The adjusted program IRR was 0.78 in the Poisson model, 0.78 in the negative binomial model, and 0.77 for the ZIP count component. The ZIP model also estimates a separate structural-zero mechanism. These are model-based associations in simulated data; a favorable fit does not prove a care-capture mechanism or establish a causal program effect.

8.1.1 Auditable reporting template

We analyzed [count definition] observed over [person-time definition] among [population and dates]. The primary estimand was [conditional rate ratio, marginal predicted count, prediction objective, or causal estimand]. We modeled the count using [Poisson/NB/ZIP] with offset(log([person-time])) and adjusted for [prespecified covariates]. We reported [IRR or other measure] with 95% confidence intervals and absolute predicted outcomes for [named covariate profiles]. We assessed [dispersion, zero frequency, functional form, dependence, calibration]. The observed zero process was interpreted as [measurement/process rationale], and limitations included [missingness, clustering, residual confounding, capture boundaries, transportability].

9 Common Errors at a Glance

Error Why it is wrong Better approach
Model a rate as ordinary Gaussian data Rates have denominator-dependent variance and can be skewed Model the count and offset its log person-time
Interpret an IRR as a risk ratio Rate and cumulative risk are distinct State the event-time scale and report absolute predictions
Treat the offset coefficient as estimated It is fixed at one by design Verify the person-time definition and use offset(log(time))
Ignore overdispersion Poisson standard errors can be too small Diagnose causes and consider quasi-Poisson, NB, clustering, or a better mean model
Choose ZIP because a histogram has many zeros Multiple mechanisms can create zeros Articulate and test the measurement/process rationale
Read ZIP μ\mu as the population mean It is conditional on the nonstructural process Report (1−π)μ(1-\pi)\mu when that is the target
Call a zero-component odds ratio an observed-zero odds ratio ZIP zero component is structural-zero membership Name the model component and scale explicitly
Compare quasi-Poisson by AIC It has no ordinary likelihood AIC Compare predictive performance or likelihood models appropriately
Ignore repeated rows or clinics Independence-based uncertainty can fail Use a design-appropriate correlated-data method
Claim causality from adjusted count regression Covariate adjustment does not solve all bias Define the causal target and identification assumptions

10 Exercises and Answers

  1. A person has two visits in 0.5 person-years. Is the outcome count, rate, or risk?
  2. Why is offset(log(person_years)) used instead of offset(person_years)?
  3. A Poisson model has Pearson dispersion 2.4. Name three possible explanations.
  4. In ZIP, what is the expected observed count for μ=1.2\mu=1.2 and π=0.25\pi=0.25?
  5. Does a lower ZIP AIC prove that some people are structurally unable to have an event?
  6. What prediction would you report for a clinician choosing between two one-year program profiles?
Show exercise answers
  1. The recorded outcome is a count of two; dividing by 0.5 gives a rate of four visits per person-year. A risk would require a specified probability of at least one event in a fixed period.
  2. The log link makes the expected count proportional to person-time: log⁡(μ)=log⁡(t)+η\log(\mu)=\log(t)+\eta. Using raw person-time as an offset imposes an implausible exponential relation in tt.
  3. Examples include omitted heterogeneity in event rates, within-clinic or within-person dependence, an incorrect functional form, an NB-like mixture, and measurement or structural-zero processes.
  4. The marginal mean is (1−0.25)×1.2=0.90(1-0.25)\times1.2=0.90 observed visits.
  5. No. AIC is a relative likelihood-based fit criterion; it cannot establish the substantive mechanism, causality, or generalizability.
  6. Report each profile’s marginal expected count and probability of at least one/zero visit, specify that follow-up is one year, and show uncertainty if the prediction will guide decisions.

11 Quick Reference

11.1 Formulas and R entry points

Target Formula or code Interpretation
Poisson rate model glm(y ~ x + offset(log(time)), family = poisson) Conditional mean count proportional to time
IRR exp⁡(βj)\exp(\beta_j) Multiplicative conditional rate comparison
Poisson zero probability e−μe^{-\mu} Probability of no event in one-process Poisson model
Pearson dispersion ∑rP2/(n−p)\sum r_P^2/(n-p) Descriptive check for extra variation
Quasi-Poisson family = quasipoisson Mean model with scale-adjusted uncertainty
Negative binomial MASS::glm.nb(...) Count model with extra-Poisson heterogeneity
ZIP zero probability π+(1−π)e−μ\pi+(1-\pi)e^{-\mu} Structural plus Poisson sources of zero
ZIP marginal mean (1−π)μ(1-\pi)\mu Expected observed count across both processes
ZIP fitting here fit_zip(count_formula, zero_formula, data) Transparent MLE with optim() and Hessian SEs

11.2 Final checklist

  • Is the event a valid nonnegative count and is the observation period explicit?
  • Is person-time positive, aligned with the count, and included as a log offset?
  • Are rate, risk, expected count, and zero probability kept distinct?
  • Are IRRs accompanied by confidence intervals and understandable absolute predictions?
  • Were dispersion, functional form, influential observations, and row dependence investigated?
  • Is the rationale for a zero mixture grounded in a measurement or event process rather than a zero bar alone?
  • Are both ZIP components and the marginal prediction scale labeled clearly?
  • Were observed-versus-predicted frequencies and calibration examined?
  • Are AIC comparisons restricted to comparable likelihood-based fits and not overinterpreted?
  • Do missing-data, clustering, selection, and causal claims match the study design?

Next directions for study

Next topics include zero-inflated negative binomial models, hurdle models, splines for nonlinear rates, distributed lags, recurrent-event survival models, generalized estimating equations, mixed-effects count models, cluster-robust inference, Bayesian hierarchical count models, and externally validated prediction models.

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.60      
##  [5] lattice_0.22-9  cachem_1.1.0    knitr_1.51      htmltools_0.5.9
##  [9] rmarkdown_2.31  stats4_4.6.1    lifecycle_1.0.5 cli_3.6.6      
## [13] grid_4.6.1      sass_0.4.10     jquerylib_0.1.4 compiler_4.6.1 
## [17] tools_4.6.1     nlme_3.1-169    evaluate_1.0.5  bslib_0.12.0   
## [21] yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0  MASS_7.3-65

For teaching only: This module introduces count-rate regression, overdispersion, zero-inflated modeling, diagnosis, and reporting. It does not replace a study-specific analysis plan, outcome-capture audit, or consultation with a biostatistician. All data and numerical results are simulated.