Zero Inflated Poisson Regression: Formula, Results, Interpretation and Diagnostics

# -*- coding: utf-8 -*-
r"""
Zero Inflated Poisson Regression - Python Output Generator

- Reads dataset.csv from the script folder first, then BASE_DIR fallback.
- Creates Python_Output/png, Python_Output/csv, Python_Output/txt, and Python_Output/pdf.
- Cleans old outputs before regenerating.
- Fits a Zero-Inflated Poisson model and a standard Poisson comparison model.
- Default count outcome: absences.
- Default count predictors: G1, G2, studytime, failures, age, health.
- Default inflation predictors: failures, age, health.
- Generates PNG charts, CSV tables, TXT report, and PDF report.
"""

from __future__ import annotations

from pathlib import Path
import os
import warnings
warnings.filterwarnings("ignore")

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.api as sm
from statsmodels.discrete.count_model import ZeroInflatedPoisson

TEST_NAME = "Zero Inflated Poisson Regression"
BASE_DIR = Path(os.environ.get(
    "ANALYSIS_BASE_DIR",
    r"analysis-folder",
))
OUTCOME_DEFAULT = "absences"
COUNT_PREDICTORS_DEFAULT = ["G1", "G2", "studytime", "failures", "age", "health"]
INFLATION_PREDICTORS_DEFAULT = ["failures", "age", "health"]
ALPHA = 0.05

OUTPUT_DIR = BASE_DIR / "Python_Output"
PNG_DIR = OUTPUT_DIR / "png"
CSV_DIR = OUTPUT_DIR / "csv"
TXT_DIR = OUTPUT_DIR / "txt"
PDF_DIR = OUTPUT_DIR / "pdf"

for folder in (OUTPUT_DIR, PNG_DIR, CSV_DIR, TXT_DIR, PDF_DIR):
    folder.mkdir(parents=True, exist_ok=True)


def cleanup_outputs() -> None:
    for folder in (PNG_DIR, CSV_DIR, TXT_DIR, PDF_DIR):
        for pattern in ("*.png", "*.csv", "*.txt", "*.pdf"):
            for item in folder.glob(pattern):
                try:
                    item.unlink()
                except OSError:
                    pass
    (OUTPUT_DIR / "00_PYTHON_SCRIPT_STARTED.txt").write_text(
        "Zero Inflated Poisson Regression Python script started successfully.\n",
        encoding="utf-8",
    )


def find_dataset(base_dir: Path) -> Path:
    script_dir = Path(__file__).resolve().parent
    candidates = [
        script_dir / "dataset.csv",
        script_dir / "dataset",
        script_dir / "dataset.txt",
        base_dir / "dataset.csv",
        base_dir / "dataset",
        base_dir / "dataset.txt",
    ]
    for candidate in candidates:
        if candidate.exists():
            return candidate
    raise FileNotFoundError(
        f"No dataset found. Expected dataset.csv in the script folder or {base_dir}."
    )


def read_dataset(path: Path) -> pd.DataFrame:
    for sep in (",", ";", "\t"):
        try:
            frame = pd.read_csv(path, sep=sep)
            if frame.shape[1] > 1:
                frame.columns = [str(c).strip() for c in frame.columns]
                return frame
        except Exception:
            continue
    frame = pd.read_csv(path)
    frame.columns = [str(c).strip() for c in frame.columns]
    return frame


def numeric_usable(series: pd.Series, min_n: int = 10) -> bool:
    values = pd.to_numeric(series, errors="coerce").dropna()
    return (
        len(values) >= min_n
        and values.nunique() > 1
        and np.isfinite(values.std())
        and values.std() > 0
    )


def choose_outcome(df: pd.DataFrame) -> str:
    if OUTCOME_DEFAULT in df.columns and numeric_usable(df[OUTCOME_DEFAULT]):
        return OUTCOME_DEFAULT
    candidates = []
    for name in df.columns:
        if not numeric_usable(df[name]):
            continue
        values = pd.to_numeric(df[name], errors="coerce").dropna()
        if np.all(values >= 0) and np.allclose(values, np.round(values)):
            candidates.append(name)
    if not candidates:
        raise ValueError("No usable non-negative count outcome was found.")
    return candidates[0]


def choose_predictors(
    df: pd.DataFrame,
    outcome: str,
    preferred: list[str],
    limit: int,
) -> list[str]:
    selected: list[str] = []
    for name in preferred:
        if name in df.columns and name != outcome and numeric_usable(df[name]):
            selected.append(name)
    for name in df.columns:
        if len(selected) >= limit:
            break
        if name == outcome or name in selected:
            continue
        if numeric_usable(df[name]):
            selected.append(name)
    if len(selected) < 1:
        raise ValueError("No usable numeric predictors were found.")
    return selected[:limit]


def prepare_data(df: pd.DataFrame):
    outcome = choose_outcome(df)
    count_predictors = choose_predictors(
        df, outcome, COUNT_PREDICTORS_DEFAULT, 6
    )
    inflation_predictors = [
        p for p in INFLATION_PREDICTORS_DEFAULT
        if p in count_predictors or (p in df.columns and numeric_usable(df[p]))
    ]
    if not inflation_predictors:
        inflation_predictors = count_predictors[: min(3, len(count_predictors))]

    columns = [outcome] + list(dict.fromkeys(count_predictors + inflation_predictors))
    data = df[columns].copy()
    for name in data.columns:
        data[name] = pd.to_numeric(data[name], errors="coerce")
    data = data.replace([np.inf, -np.inf], np.nan).dropna().reset_index(drop=True)

    data[outcome] = np.round(data[outcome]).astype(int)
    data = data[data[outcome] >= 0].reset_index(drop=True)
    if len(data) <= len(count_predictors) + len(inflation_predictors) + 10:
        raise ValueError("Not enough complete observations for Zero-Inflated Poisson Regression.")
    if int((data[outcome] == 0).sum()) == 0:
        raise ValueError("The selected count outcome contains no zero observations.")

    return data, outcome, count_predictors, inflation_predictors


def fit_models(data, outcome, count_predictors, inflation_predictors):
    y = data[outcome].astype(int)
    x_count = sm.add_constant(
        data[count_predictors].astype(float), has_constant="add"
    )
    x_infl = sm.add_constant(
        data[inflation_predictors].astype(float), has_constant="add"
    )

    poisson = sm.GLM(y, x_count, family=sm.families.Poisson()).fit()

    zip_model = ZeroInflatedPoisson(
        endog=y,
        exog=x_count,
        exog_infl=x_infl,
        inflation="logit",
    )
    methods = ["bfgs", "lbfgs", "newton"]
    zip_result = None
    last_error = None
    for method in methods:
        try:
            candidate = zip_model.fit(
                method=method,
                maxiter=1500,
                disp=0,
            )
            zip_result = candidate
            converged = bool(candidate.mle_retvals.get("converged", True))
            if converged:
                break
        except Exception as exc:
            last_error = exc
    if zip_result is None:
        raise RuntimeError(f"Zero-Inflated Poisson fitting failed: {last_error}")

    return y, x_count, x_infl, poisson, zip_result


def split_coefficients(zip_result, count_predictors, inflation_predictors):
    params = zip_result.params
    bse = zip_result.bse
    pvalues = zip_result.pvalues
    conf = zip_result.conf_int(alpha=ALPHA)

    count_names = ["const"] + count_predictors
    inflation_names = ["inflate_const"] + [f"inflate_{p}" for p in inflation_predictors]

    count_rows = []
    for name in count_names:
        beta = float(params[name])
        se = float(bse[name])
        p = float(pvalues[name])
        low = float(conf.loc[name, 0])
        high = float(conf.loc[name, 1])
        count_rows.append({
            "term": "Intercept" if name == "const" else name,
            "coefficient_log_count": beta,
            "standard_error": se,
            "z_value": beta / se,
            "p_value": p,
            "ci95_lower_log": low,
            "ci95_upper_log": high,
            "incidence_rate_ratio": np.exp(beta),
            "irr_ci95_lower": np.exp(low),
            "irr_ci95_upper": np.exp(high),
            "decision_alpha_0_05": (
                "Statistically significant"
                if p < ALPHA else "Not statistically significant"
            ),
        })

    inflation_rows = []
    for name in inflation_names:
        beta = float(params[name])
        se = float(bse[name])
        p = float(pvalues[name])
        low = float(conf.loc[name, 0])
        high = float(conf.loc[name, 1])
        display = "Intercept" if name == "inflate_const" else name.replace("inflate_", "")
        inflation_rows.append({
            "term": display,
            "coefficient_logit_zero": beta,
            "standard_error": se,
            "z_value": beta / se,
            "p_value": p,
            "ci95_lower_logit": low,
            "ci95_upper_logit": high,
            "odds_ratio_structural_zero": np.exp(beta),
            "or_ci95_lower": np.exp(low),
            "or_ci95_upper": np.exp(high),
            "decision_alpha_0_05": (
                "Statistically significant"
                if p < ALPHA else "Not statistically significant"
            ),
        })

# Zero Inflated Poisson Regression - R Output Generator
#
# - Reads dataset.csv from the script folder first, then BASE_DIR fallback.
# - Creates R_Output/charts, R_Output/tables, R_Output/txt, and R_Output/pdf.
# - Cleans old outputs before regenerating.
# - Fits a Zero-Inflated Poisson model using pscl::zeroinfl.
# - Fits a standard Poisson comparison model.
# - Generates colorful PNG charts, CSV tables, TXT report, and PDF report.

BASE_DIR <- Sys.getenv(
  "ANALYSIS_BASE_DIR",
  unset = "analysis-folder"
)

TEST_NAME <- "Zero Inflated Poisson Regression"
OUTCOME_DEFAULT <- "absences"
COUNT_PREDICTORS_DEFAULT <- c("G1", "G2", "studytime", "failures", "age", "health")
INFLATION_PREDICTORS_DEFAULT <- c("failures", "age", "health")
ALPHA <- 0.05

OUTPUT_DIR <- file.path(BASE_DIR, "R_Output")
CHART_DIR <- file.path(OUTPUT_DIR, "charts")
TABLE_DIR <- file.path(OUTPUT_DIR, "tables")
TXT_DIR <- file.path(OUTPUT_DIR, "txt")
PDF_DIR <- file.path(OUTPUT_DIR, "pdf")

for (folder in c(OUTPUT_DIR, CHART_DIR, TABLE_DIR, TXT_DIR, PDF_DIR)) {
  if (!dir.exists(folder)) dir.create(folder, recursive = TRUE, showWarnings = FALSE)
}

COLORS <- c(
  "#2E86AB", "#F18F01", "#C73E1D", "#6A994E",
  "#6F4E7C", "#4D908E", "#F9C74F", "#577590"
)

cleanup_outputs <- function() {
  for (folder in c(CHART_DIR, TABLE_DIR, TXT_DIR, PDF_DIR)) {
    for (pattern in c("*.png", "*.csv", "*.txt", "*.pdf")) {
      files <- Sys.glob(file.path(folder, pattern))
      if (length(files) > 0) unlink(files, force = TRUE)
    }
  }
  writeLines(
    "Zero Inflated Poisson Regression R script started successfully.",
    file.path(OUTPUT_DIR, "00_R_SCRIPT_STARTED.txt")
  )
}

find_dataset <- function(base_dir) {
  candidates <- c()
  args <- commandArgs(trailingOnly = FALSE)
  match <- grep("^--file=", args)
  if (length(match) > 0) {
    script_path <- normalizePath(
      sub("^--file=", "", args[match[1]]),
      winslash = "/",
      mustWork = FALSE
    )
    script_dir <- dirname(script_path)
    candidates <- c(
      candidates,
      file.path(script_dir, "dataset.csv"),
      file.path(script_dir, "dataset"),
      file.path(script_dir, "dataset.txt")
    )
  }
  candidates <- c(
    candidates,
    file.path(base_dir, "dataset.csv"),
    file.path(base_dir, "dataset"),
    file.path(base_dir, "dataset.txt")
  )
  for (item in candidates) if (file.exists(item)) return(item)
  stop("dataset.csv was not found in the script folder or BASE_DIR.")
}

read_dataset <- function(path) {
  for (sep in c(",", ";", "\t")) {
    attempt <- try(
      read.csv(path, sep = sep, stringsAsFactors = FALSE, check.names = FALSE),
      silent = TRUE
    )
    if (!inherits(attempt, "try-error") && ncol(attempt) > 1) {
      names(attempt) <- trimws(names(attempt))
      return(attempt)
    }
  }
  out <- read.csv(path, stringsAsFactors = FALSE, check.names = FALSE)
  names(out) <- trimws(names(out))
  out
}

numeric_usable <- function(x, min_n = 10) {
  values <- suppressWarnings(as.numeric(x))
  values <- values[is.finite(values)]
  length(values) >= min_n &&
    length(unique(values)) > 1 &&
    is.finite(sd(values)) &&
    sd(values) > 0
}

prepare_data <- function(df) {
  outcome <- if (
    OUTCOME_DEFAULT %in% names(df) &&
    numeric_usable(df[[OUTCOME_DEFAULT]])
  ) OUTCOME_DEFAULT else {
    candidates <- names(df)[sapply(df, function(x) {
      if (!numeric_usable(x)) return(FALSE)
      v <- suppressWarnings(as.numeric(x))
      v <- v[is.finite(v)]
      all(v >= 0) && all(abs(v - round(v)) < 1e-8)
    })]
    if (length(candidates) == 0) stop("No usable count outcome was found.")
    candidates[1]
  }

  count_predictors <- c()
  for (name in COUNT_PREDICTORS_DEFAULT) {
    if (name %in% names(df) && name != outcome && numeric_usable(df[[name]])) {
      count_predictors <- c(count_predictors, name)
    }
  }
  for (name in names(df)) {
    if (length(count_predictors) >= 6) break
    if (name == outcome || name %in% count_predictors) next
    if (numeric_usable(df[[name]])) count_predictors <- c(count_predictors, name)
  }

  inflation_predictors <- INFLATION_PREDICTORS_DEFAULT[
    INFLATION_PREDICTORS_DEFAULT %in% names(df)
  ]
  inflation_predictors <- inflation_predictors[
    sapply(inflation_predictors, function(name) numeric_usable(df[[name]]))
  ]
  if (length(inflation_predictors) == 0) {
    inflation_predictors <- head(count_predictors, min(3, length(count_predictors)))
  }

  columns <- unique(c(outcome, count_predictors, inflation_predictors))
  data <- df[, columns, drop = FALSE]
  for (name in names(data)) data[[name]] <- suppressWarnings(as.numeric(data[[name]]))
  data <- data[complete.cases(data), , drop = FALSE]
  data[[outcome]] <- round(data[[outcome]])
  data <- data[data[[outcome]] >= 0, , drop = FALSE]
  rownames(data) <- NULL

  if (sum(data[[outcome]] == 0) == 0) {
    stop("The selected outcome contains no zero observations.")
  }

  list(
    data = data,
    outcome = outcome,
    count_predictors = unique(count_predictors),
    inflation_predictors = unique(inflation_predictors)
  )
}

safe_name <- function(x) paste0("`", gsub("`", "", x), "`")

make_formula <- function(outcome, count_predictors, inflation_predictors) {
  count_rhs <- paste(safe_name(count_predictors), collapse = " + ")
  infl_rhs <- paste(safe_name(inflation_predictors), collapse = " + ")
  as.formula(paste(
    safe_name(outcome), "~", count_rhs, "|", infl_rhs
  ))
}

coefficient_tables <- function(model) {
  s <- summary(model)
  count <- as.data.frame(s$coefficients$count)
  count$term <- rownames(count)
  rownames(count) <- NULL
  names(count)[1:4] <- c("coefficient_log_count", "standard_error", "z_value", "p_value")
  count$ci95_lower_log <- count$coefficient_log_count - qnorm(0.975) * count$standard_error
  count$ci95_upper_log <- count$coefficient_log_count + qnorm(0.975) * count$standard_error
  count$incidence_rate_ratio <- exp(count$coefficient_log_count)
  count$irr_ci95_lower <- exp(count$ci95_lower_log)
  count$irr_ci95_upper <- exp(count$ci95_upper_log)
  count$decision_alpha_0_05 <- ifelse(
    count$p_value < ALPHA,
    "Statistically significant",
    "Not statistically significant"
  )
  count$term <- ifelse(count$term == "(Intercept)", "Intercept", gsub("`", "", count$term))
  count <- count[, c(
    "term", "coefficient_log_count", "standard_error", "z_value", "p_value",
    "ci95_lower_log", "ci95_upper_log", "incidence_rate_ratio",
    "irr_ci95_lower", "irr_ci95_upper", "decision_alpha_0_05"
  )]

  infl <- as.data.frame(s$coefficients$zero)
  infl$term <- rownames(infl)
  rownames(infl) <- NULL
  names(infl)[1:4] <- c("coefficient_logit_zero", "standard_error", "z_value", "p_value")
  infl$ci95_lower_logit <- infl$coefficient_logit_zero - qnorm(0.975) * infl$standard_error
  infl$ci95_upper_logit <- infl$coefficient_logit_zero + qnorm(0.975) * infl$standard_error
  infl$odds_ratio_structural_zero <- exp(infl$coefficient_logit_zero)
  infl$or_ci95_lower <- exp(infl$ci95_lower_logit)
  infl$or_ci95_upper <- exp(infl$ci95_upper_logit)
  infl$decision_alpha_0_05 <- ifelse(
    infl$p_value < ALPHA,
    "Statistically significant",
    "Not statistically significant"
  )
  infl$term <- ifelse(infl$term == "(Intercept)", "Intercept", gsub("`", "", infl$term))
  infl <- infl[, c(
    "term", "coefficient_logit_zero", "standard_error", "z_value", "p_value",
    "ci95_lower_logit", "ci95_upper_logit", "odds_ratio_structural_zero",
    "or_ci95_lower", "or_ci95_upper", "decision_alpha_0_05"
  )]

  list(count = count, inflation = infl)
}

expected_frequencies <- function(model, poisson, data, outcome) {
  y <- data[[outcome]]
  max_count <- max(y)
  lambda <- predict(model, type = "count")
  pi_zero <- predict(model, type = "zero")
  poisson_mu <- fitted(poisson)

  rows <- lapply(0:max_count, function(k) {
    zip_prob <- if (k == 0) {
      pi_zero + (1 - pi_zero) * dpois(0, lambda)
    } else {
      (1 - pi_zero) * dpois(k, lambda)
    }
    data.frame(
      count_value = k,
      observed_frequency = sum(y == k),
      expected_zip_frequency = sum(zip_prob),
      expected_poisson_frequency = sum(dpois(k, poisson_mu)),
      stringsAsFactors = FALSE
    )
  })
  do.call(rbind, rows)
}

save_chart <- function(filename, plot_function, width = 11, height = 7) {
  path <- file.path(CHART_DIR, filename)
  png(path, width = width, height = height, units = "in", res = 180)
  tryCatch(
    plot_function(),
    error = function(e) {
      plot.new()
      title(main = paste("Chart generation failed:", filename))
      text(0.05, 0.85, e$message, adj = c(0, 1))
    },
    finally = dev.off()
  )
  path
}

create_charts <- function(
  data, outcome, count_coef, infl_coef, comparison, frequency,
  predictions, model_summary
) {
  paths <- c()

  paths <- c(paths, save_chart("01_observed_count_distribution.png", function() {
    counts <- table(data[[outcome]])
    barplot(
      counts, col = COLORS[1], border = "#1F2937",
      xlab = outcome, ylab = "Observed frequency",

Standard Poisson assumes one count-generating process and assigns every zero to the Poisson distribution. Zero Inflated Poisson Regression adds a second state that always produces zero. In the worked data, Poisson expects only 24.6591 zeros, whereas ZIP expects 243.8578 against 244 observed.

The AIC difference of 1286.2765 strongly favors ZIP for the fitted distribution. The comparison does not show that every aspect of ZIP is adequate; it shows that the one-process Poisson model is unable to represent the zero mass.

Negative binomial regression handles overdispersion by adding a dispersion parameter to one count process. It can produce more zeros than Poisson indirectly because its variance is larger, but it does not posit a separate always-zero state.

Compare Zero Inflated Poisson Regression with negative binomial regression when the substantive case for structural zeros is uncertain. A negative binomial model may be more parsimonious if extra zeros are simply one consequence of broad count heterogeneity.

Zero-inflated negative binomial regression retains the structural-zero mixture and replaces the Poisson count state with a negative binomial state. It is designed for data containing both excess zeros and extra-Poisson variation among positive counts.

The variance-to-mean ratio of 5.8852 and the long positive residual tail make ZINB a particularly important sensitivity model here. Zero Inflated Poisson Regression should not be selected as final until the dispersion parameter and information criteria from ZINB are examined.

A hurdle model uses a binary equation for zero versus positive and a zero-truncated count model for positive outcomes. Therefore, every zero belongs to the first process. Zero Inflated Poisson Regression allows zeros from both the structural state and the Poisson count state.

The choice is substantive as well as statistical. Use a hurdle model when crossing from zero to any positive count is a distinct event and positive counts follow a separate process. Use ZIP when count-state cases can legitimately produce zero.

A structural zero is generated by the additional zero state. A sampling zero is generated by the Poisson state when the random count happens to equal zero. The same observed value can arise through either route.

Zero Inflated Poisson Regression estimates probabilities, not observed labels. Statements such as “244 structural zeros” are incorrect. The verified result is 244 observed zeros and 243.8578 expected total zeros under the fitted mixture.

Inflation predictors should explain why some cases have a qualitatively different chance of always producing zero. They do not need to match the count predictors, but the choice should follow theory, design and availability before model fitting.

In this analysis, failures, age and health enter the inflation equation. Their individual p-values exceed .05, and the predicted zero probabilities overlap strongly across outcome groups. Alternative predictors may be needed to explain latent zero-state membership.

Count predictors explain variation in λ among cases in the count-generating state. They can include continuous variables, dummy-coded categories, interactions, splines and an exposure offset.

The worked count equation includes G1, G2, studytime, failures, age and health. Zero Inflated Poisson Regression interprets their exponentiated coefficients as IRRs conditional on count-state membership.

The count intercept is the log expected count when every count predictor equals zero. Exponentiating 0.6639 gives 1.9423, but the corresponding profile is unrealistic because age and grade variables do not naturally equal zero in the analysed range.

Centering continuous predictors at meaningful values makes the intercept useful without changing slopes or fitted probabilities. For example, centering age and grades at their sample means would define the baseline λ for an average profile.

The inflation intercept is the log odds of structural-zero membership when all inflation predictors equal zero. Its OR is 2.7109, but the 95% interval is wide and includes one.

As with the count intercept, raw zero values may not define a meaningful student profile. Centering age and health can improve interpretation and numerical stability while preserving the Zero Inflated Poisson Regression fit.

The inflation probability π is only the chance of the always-zero state. The complete probability of observing zero adds the chance of entering the count state and drawing a Poisson zero.

Reporting π as P(Y = 0) understates zero probability. The charts and calibration table in this post use the full mixture probability, which has an all-case mean of 0.3757.

The final fitted count is μ = (1 − π)λ. It is lower than λ whenever the structural-zero probability is positive. This quantity should be used when plotting observed counts against fitted values.

Zero Inflated Poisson Regression software often returns several prediction types. Analysts must document whether a column contains η, λ, π, μ or total zero probability before calculating residuals or accuracy metrics.

A count IRR describes λ, and an inflation OR describes π. Neither is automatically the effect on μ. Observed-scale marginal effects combine both equations and can vary across predictor profiles.

For a predictor appearing in both equations, the count and inflation pathways may reinforce or offset one another. Average marginal effects should therefore be calculated directly rather than inferred from one transformed coefficient.

Count models become rate models when cases have different observation times or opportunities. The standard remedy is an offset equal to the log of exposure in the count equation with its coefficient fixed at one.

The worked Zero Inflated Poisson Regression contains no offset and assumes comparable absence opportunities. Applying the model to unequal follow-up periods without exposure adjustment would confound event intensity with time at risk.

Zero inflation can reduce apparent overdispersion because it adds probability at zero, but the count state still obeys the Poisson mean-variance equality. Residual overdispersion can remain among positive counts.

The large positive Pearson residuals and underprediction of high absences show that Zero Inflated Poisson Regression has not completely solved dispersion in this dataset. ZINB is the natural next model.

The observed frequency distribution has large spikes at several even counts and much smaller frequencies at adjacent odd counts. Such heaping may reflect reporting, scheduling or rounding rather than the underlying event process.

A smooth Poisson mixture cannot reproduce a heaping mechanism exactly. Analysts should inspect data collection rules and consider grouped-count or heaping models before attributing every mismatch to omitted predictors.

ZIP estimation must distinguish the inflation state from ordinary low Poisson means. Identification is stronger when predictors explain the two components differently and when the data contain informative variation in zeros and positive counts.

The narrow range of predicted zero probabilities suggests weak separation in this example. Convergence alone does not prove that the latent classes have a clear substantive interpretation.

Binary logit components can suffer from separation when a predictor or combination perfectly predicts zero-state membership. Estimates may become extremely large, standard errors may inflate and optimization may become unstable.

No such extreme coefficients appear in the supplied Zero Inflated Poisson Regression, but the inflation equation should still be checked for sparse categories, extreme fitted probabilities and sensitivity to predictor selection.

Highly correlated predictors can inflate standard errors and make coefficient signs unstable in either component. G1 and G2 are conceptually related, so their conditional interpretations depend on including both.

Variance-inflation diagnostics can be computed separately for the count and inflation design matrices. Removing a variable solely because it is nonsignificant can change the estimand and should be guided by theory rather than an automatic threshold.

The log and logit equations assume linear predictor effects on their respective link scales. Age, grade and health relationships may be nonlinear, and one predictor may modify another.

Splines, polynomial terms and prespecified interactions can improve Zero Inflated Poisson Regression when supported by subject knowledge and sufficient data. Comparisons should use likelihood, information criteria, residual plots and validation rather than p-values alone.

Students nested within schools or measured repeatedly are not independent. Unmodelled dependence can distort standard errors, produce extra variation and mimic zero inflation.

Cluster-robust covariance estimates, generalized estimating equations, random-effects count models or multilevel zero-inflated models may be required. The present analysis treats the 649 records as independent.

Pearson residuals standardize observed-minus-predicted differences using model-based variance. Randomized quantile residuals can provide a more continuous diagnostic for discrete mixture models.

For this Zero Inflated Poisson Regression, the long positive residual tail is more important than minor central scatter. It points to high-count underfit and should trigger case review and alternative variance models.

Calibration asks whether predicted probabilities agree with observed frequencies. Discrimination asks whether higher probabilities are assigned to cases that actually experience the outcome.

The model is strongly calibrated for total zeros but weakly discriminative: mean zero probability is 0.3809 for observed zeros and 0.3726 for positive counts. Both findings should be reported.

AIC and BIC compare the tradeoff between fit and complexity. The Vuong statistic compares average case-level log-likelihood contributions for nonnested candidates under specific assumptions.

ZIP strongly outperforms Poisson on all three measures, but these comparisons do not certify correct specification. They should be combined with residual diagnostics and comparisons involving negative-binomial alternatives.

In-sample likelihood and error metrics are optimistic for new data. Proper validation requires refitting on training data and evaluating μ, P(Y = 0), positive-count error and calibration on untouched observations.

Because Zero Inflated Poisson Regression has two prediction targets, validation should not collapse everything into one RMSE. Zero calibration and positive-count forecasting answer different practical questions.

Complete-case analysis is valid only under assumptions about why values are missing. Different missing-data handling can change both the zero proportion and the relationships in each component.

A reproducible report should state the original sample size, exclusions, formulas, software versions, optimizer, convergence status, coefficient parameterization and prediction type. The downloadable reports and workbook preserve the fitted calculations for this example.

Case-level output should place observed counts beside λ, π, μ, total zero probability and residuals. This makes it possible to identify high-count underprediction, unusual zero probabilities and records that dominate fit.

The worked Excel file provides the practical audit route without placing hundreds of rows inside the WordPress article. Review extreme Pearson residuals, compare fitted components, and document any data correction or sensitivity refit rather than silently deleting unusual observations.