Package {underdisp}

underdisp: Diagnostics and Models for Underdispersed Count Data

Description

Detect and model underdispersion (conditional variance below the conditional mean) in count data. The package provides a screening diagnostic (ud_screen()), the continuous parameter binomial regression model (cpb()) with a zero-truncated variant, validated bootstrap and profile-likelihood inference, a family-comparison helper (compare_dispersion()), and King-style quantities of interest (predict.cpb(), implied_ceiling(), first_difference()).

Author(s)

Maintainer: Benjamin E. Bagozzi bagozzib@udel.edu

See Also

Useful links:

• https://github.com/bagozzib/underdisp

• Report bugs at https://github.com/bagozzib/underdisp/issues

Profile-likelihood interval for the dispersion parameter alpha

Description

Profile-likelihood interval for the dispersion parameter alpha

Usage

alpha_confint(object, level = 0.95)

Arguments

Value

A length-2 numeric vector (lower, upper) with attributes alpha (the point estimate) and boundary (TRUE if the lower bound is at the feasibility boundary, i.e. strong underdispersion, where the interval is one-sided).

Examples

set.seed(7); x <- rnorm(300)
N <- pmax(round(exp(1.5 + 0.4 * x) / 0.5), 1); y <- rbinom(300, N, 0.5)
fit <- cpb(y ~ x, data.frame(y = y, x = x)[y > 0, ], se = "none")
alpha_confint(fit)

Augment data with CPB fitted values (broom method)

Description

Augment data with CPB fitted values (broom method)

Usage

## S3 method for class 'cpb'
augment(x, ...)

Arguments

Value

A data frame with fitted values, response residuals, and the implied per-observation ceiling.

Compare count-model dispersion across the model family

Description

Fits the Poisson and negative binomial defaults alongside the two underdispersed workhorses—the soft-tail Conway–Maxwell–Poisson (fit natively, no external dependency) and the hard-ceiling CPB—and returns a comparison on both fit (log-likelihood, AIC, BIC) and calibration (mean logarithmic score and ranked probability score, lower is better, plus the fitted share of zeros against the observed share). Whether the tail is better described by an accelerated decay (COM-Poisson) or a binding ceiling (CPB) is a testable question, not an assumption; the scoring rules are the calibration counterpart to the information criteria. The CPB's ceiling-exceedance share (observations above the implied ceiling) is reported as a falsification check on the hard bound.

Usage

compare_dispersion(
  formula,
  data,
  max.support = 500,
  hurdle = FALSE,
  zi = FALSE
)

Arguments

Details

With hurdle = TRUE and zeros present, a hurdle-CPB is added, so a zero-inflated underdispersed process can be compared to the single-equation models on the same footing.

Value

A list with table (a data frame of df, logLik, AIC, BIC, logscore, rps, and zero_fit per model), obs_zero (the observed zero share), nu (the COM-Poisson dispersion, ⁠>1⁠ = underdispersion, NA if the COM-Poisson fit failed), alpha (the CPB shape parameter), ceiling_exceedance (share of observations above the CPB ceiling), and cpb_ok (FALSE if the single-equation CPB could not satisfy its feasibility constraint on the data, e.g. under heavy zero-inflation with a wide count range — itself a signal that a hurdle or zero-inflated model is needed).

See Also

ud_screen(), cpb(), score()

Examples

set.seed(1); x <- rnorm(120)
N <- pmax(round(exp(1.6 + 0.4 * x) / 0.5), 1); y <- rbinom(120, N, 0.5)
compare_dispersion(y ~ x, data = data.frame(y = y, x = x))$table

Compare fitted underdispersed-count models

Description

Compares fitted models from this package — cpb, cpb_fe, hurdle_cpb, or zi_cpb, in any combination of fixed effects and robust/clustered standard errors — on information criteria and, where a predicted distribution is available, on proper scores. Unlike zi_test() this is not a hypothesis test: the models need not be nested, so it is the appropriate tool for the non-nested hurdle-versus-mixture choice. Passing an object that is not a fitted model from this package is an error, and models fit on different data raise a warning.

Usage

compare_models(...)

Arguments

Value

A data frame with df, logLik, AIC, BIC, and (where the predicted distribution is available) logscore and rps, one row per model, ordered by AIC.

See Also

zi_test(), compare_dispersion()

Examples

set.seed(1); n <- 800; x <- rnorm(n); z <- rnorm(n)
y <- rhurdle_cpb(n, exp(1.2 + 0.5 * x), 0.5, plogis(-0.2 + 0.8 * z))
d <- data.frame(y = y, x = x, z = z)
compare_models(hurdle = hurdle_cpb(y ~ x, data = d, participation = ~ z),
               zi = zi_cpb(y ~ x, data = d, zero = ~ z))

Conway–Maxwell–Poisson distribution functions

Description

Density, distribution, quantile, and random generation for the COM-Poisson with rate lambda and dispersion nu (nu > 1 underdispersed, nu = 1 Poisson, nu < 1 overdispersed). Complements the estimator count_reg(..., family = "compois").

Usage

dcompois(x, lambda, nu, log = FALSE)

pcompois(q, lambda, nu, lower.tail = TRUE, log.p = FALSE)

qcompois(p, lambda, nu, lower.tail = TRUE, log.p = FALSE)

rcompois(n, lambda, nu)

Arguments

Value

dcompois a density, pcompois a CDF, qcompois a quantile, rcompois a numeric vector of count draws.

Examples

dcompois(0:5, lambda = 3, nu = 1.5)
mean(rcompois(1000, lambda = 3, nu = 1.5))

Confidence intervals for a CPB fit

Description

Coefficient intervals use the cold-multistart bootstrap percentile method (validated to nominal coverage); the interval for alpha uses the profile-likelihood method, which is reliable except under strong underdispersion, where alpha sits at the feasibility boundary and the interval is one-sided.

Usage

## S3 method for class 'cpb'
confint(object, parm, level = 0.95, ...)

Arguments

Value

A matrix of lower/upper bounds.

Poisson and negative-binomial count regression (matched CPB baselines)

Description

Fits a Poisson or negative-binomial regression with a log link, optionally zero-truncated and/or with unit fixed effects, returning an object that compare_models(), score(), and the broom methods treat on the same footing as a cpb() or gec() fit. This is the matched baseline for the underdispersed models: it exists so a Poisson/NB can be compared to a CPB with identical degrees-of-freedom, log-likelihood, proper-score, and robust-standard-error accounting, rather than reconciled across packages.

Usage

count_reg(
  formula,
  data,
  family = c("poisson", "negbin", "compois"),
  truncated = FALSE,
  fe = NULL,
  offset = NULL,
  se = c("analytic", "robust", "cluster", "none"),
  cluster = NULL,
  ...
)

Arguments

Value

An object of class "count_reg". The component ⁠$theta⁠ holds the shape parameter: the negative-binomial size for "negbin", the COM-Poisson dispersion \nu for "compois", and NA for "poisson".

See Also

cpb(), compare_models(), hurdle_count(), zi_count()

Examples

set.seed(1); n <- 400; x <- rnorm(n)
y <- rpois(n, exp(1 + 0.5 * x))
m <- count_reg(y ~ x, data = data.frame(y = y, x = x), family = "poisson")
compare_models(pois = m,
               nb = count_reg(y ~ x, data = data.frame(y = y, x = x), family = "negbin"))

Fit a Continuous Parameter Binomial (CPB) regression

Description

Fits the underdispersed continuous parameter binomial model of King (1989), in which the conditional variance is a fraction of the conditional mean, \mathrm{Var}(Y\mid x) = \alpha\,\mathrm{E}(Y\mid x) with 0 < \alpha < 1, and each observation has an endogenous ceiling \lambda_i/(1-\alpha). A zero-truncated variant (the default) conditions on Y \ge 1, appropriate when underdispersion lives among the positive counts of an otherwise zero-inflated outcome.

Usage

cpb(
  formula,
  data,
  truncated = TRUE,
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  offset = NULL,
  alpha.start = 0.5,
  max.support = 500,
  maxit = 20000,
  reltol = 1e-08
)

Arguments

Details

The mean is modelled log-linearly, \lambda_i = \exp(x_i'\beta). Because the support depends on the parameters, the log-likelihood is non-smooth at the feasibility boundary and the numerical Hessian is unreliable; inference therefore uses a cold-multistart bootstrap for the coefficients (validated to nominal coverage) and a profile-likelihood interval for \alpha (see confint.cpb()).

Value

An object of class "cpb": a list with coefficients, alpha, bootstrap se.beta/vcov/ci.beta, loglik, loglik.null, fitted.values, ceiling (the observation-specific implied ceiling), the model frame pieces, and the original call.

References

King, G. (1989). Variance specification in event count models. American Journal of Political Science, 33(3), 762-784.

See Also

ud_screen(), confint.cpb(), predict.cpb(), implied_ceiling()

Examples

set.seed(1)
n <- 300; x <- rnorm(n)
N <- pmax(round(exp(1.6 + 0.5 * x) / 0.5), 1)
y <- rbinom(n, N, 0.5)                 # underdispersed (var/mean approx 0.5)
d <- data.frame(y = y, x = x)
fit <- cpb(y ~ x, data = d[d$y > 0, ], se = "none")
summary(fit)

Continuous parameter binomial distribution functions

Description

Density, distribution function, and quantile function for the continuous parameter binomial (CPB), and its zero-truncated form. The CPB has a hard ceiling at \lfloor \lambda/(1-\alpha) \rfloor; probability above it is zero. Complements the simulator rcpb().

Usage

dcpb(x, lambda, alpha, truncated = FALSE, log = FALSE)

pcpb(q, lambda, alpha, truncated = FALSE, lower.tail = TRUE, log.p = FALSE)

qcpb(p, lambda, alpha, truncated = FALSE, lower.tail = TRUE, log.p = FALSE)

Arguments

Value

dcpb a density, pcpb a distribution function, qcpb a quantile.

See Also

rcpb()

Examples

dcpb(0:5, lambda = 3, alpha = 0.5)
pcpb(3, lambda = 3, alpha = 0.5)
qcpb(0.9, lambda = 3, alpha = 0.5)

CPB regression with high-dimensional unit fixed effects

Description

Fits the CPB with a full set of unit fixed effects by concentrating (profiling) out the unit intercepts. Each unit's intercept is solved by a one-dimensional inner maximization, so the outer optimizer handles only the covariate coefficients and the dispersion parameter. This avoids an explicit unit-dummy design matrix and scales to thousands of units, where dummy-based fitting exhausts memory. The concentrated log-likelihood equals the full-dummy log-likelihood exactly.

Usage

cpb_fe(
  formula,
  data,
  fe,
  truncated = FALSE,
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  offset = NULL,
  max.support = 500L,
  inner_it = 30L,
  maxit = 3000L,
  reltol = 1e-07,
  bias_correct = c("none", "jackknife")
)

Arguments

Details

Short panels: the dispersion parameter alpha and the fixed effects are subject to the incidental-parameters bias of nonlinear fixed-effects estimation. The bias in alpha is downward and of order 1/T, where T is the per-unit number of observations. It is small for T at least about 30 and should be treated cautiously for short panels. The covariate coefficients are not materially affected. bias_correct = "jackknife" removes the leading 1/T term by the split-panel (half-panel) jackknife: the model is refit on the first and second temporal halves of every unit's series (rows are split in the order supplied, which should be temporal order) and the corrected estimate is 2 * full - mean(halves). In the package's fixed-effects bias Monte Carlo the alpha bias at T = 6/10/20/40 falls from -0.147/-0.096/-0.056/-0.033 to -0.038/-0.022/-0.016/-0.009. With correction on, coefficients and alpha are the corrected estimates (the maximum-likelihood values are kept in ⁠$uncorrected⁠), the unit effects and fitted values are re-concentrated at the corrected parameters, and logLik/AIC continue to refer to the maximum-likelihood fit.

Validity gate: the split-panel identity requires the two half-panels to estimate the same pseudo-true parameter (time-homogeneity; Dhaene & Jochmans 2015, Review of Economic Studies). On trending or time-heterogeneous panels the correction is invalid, so the function REFUSES it – returning the uncorrected maximum-likelihood fit with a warning naming the failed check – whenever the corrected dispersion leaves its feasible space or the two half-panel dispersion estimates disagree beyond what the 1/T bias can explain. A refusal is diagnostic information about the panel, not an error. The gate is deliberately powered over sized: in the package's calibration Monte Carlo it refuses about 9 percent of genuinely time-homogeneous panels (a conservative nuisance; the returned fit is exactly the ordinary maximum-likelihood estimate) while catching 98 percent of dispersion regime changes and 100 percent of smooth unmodeled trends – the cases where an uncaught correction would be silently wrong.

One set of fixed effects is concentrated out. For two-way (e.g. unit and time) fixed effects, add factor(time) to formula (entered as dummies), or use count_reg(), which supports two-way fixed effects natively.

Value

An object of class "cpb_fe" with coefficients, alpha, loglik, fe (the estimated unit intercepts), fitted.values, and the implied per-observation ceiling.

See Also

cpb(), ud_screen()

Examples

set.seed(1)
d <- do.call(rbind, lapply(1:60, function(i) {
  x <- rnorm(20); lam <- exp(rnorm(1, 0, 0.5) + 0.5 * x)
  N <- pmax(round(lam / 0.5), 1)
  data.frame(unit = i, x = x, y = rbinom(20, N, 0.5))
}))
fit <- cpb_fe(y ~ x, data = d, fe = "unit")
fit

Cross-validated proper scores for a count model

Description

K-fold cross-validated log score and RPS: the data are split into k folds, the model is refit on each training set via fitfun, and its predicted distribution is scored on the held-out fold. Because the score is genuinely out of sample it penalizes over-parameterization, so it is the honest criterion for comparing models that differ in the number of parameters (for example a hurdle with per-unit participation fixed effects against a zero-inflated model). Per-observation held-out scores are returned so two models fit on the same folds can be compared with a paired, cluster-robust test.

Usage

cv_score(fitfun, data, k = 5, kmax = NULL, folds = NULL)

Arguments

Value

An object of class "cv_score": a list with the mean held-out logscore and rps, the per-observation vectors logscore_i/rps_i, and the folds used.

See Also

score()

Examples

set.seed(1)
d <- data.frame(y = rcpb(500, lambda = 3, alpha = 0.5))
cv <- cv_score(function(tr) cpb(y ~ 1, tr, truncated = FALSE, se = "none"), d, k = 5)
cv$logscore

King-style first difference for a CPB fit

Description

The effect on a quantity of interest of moving one covariate from one value to another, holding the other covariates at their means, with a bootstrap percentile interval (Tomz, Wittenberg & King style, using the model's bootstrap draws).

Usage

first_difference(object, ...)

## S3 method for class 'cpb'
first_difference(
  object,
  variable,
  from,
  to,
  quantity = c("mean", "ceiling", "prob"),
  y = NULL,
  level = 0.95,
  ...
)

Arguments

Value

A "ud_fd" data frame – the package-wide first-difference contract (columns component, from, to, diff, lower, upper, method) – with one row for the requested quantity and a bootstrap percentile interval on the difference (method = "bootstrap (stored)"). Every first_difference() method in the package returns this same shape.

Examples

set.seed(1); x <- rnorm(400)
N <- pmax(round(exp(1.6 + 0.5 * x) / 0.5), 1); y <- rbinom(400, N, 0.5)
fit <- cpb(y ~ x, data = data.frame(y = y, x = x)[y > 0, ], se = "bootstrap", B = 200)
first_difference(fit, "x", from = -1, to = 1, quantity = "mean")

First differences for the matched count families

Description

The discrete-change effect on the expected count E(Y) as variable moves from from to to, holding the other covariates at their sample means. For count_reg this is a single number with a delta-method interval; for the two-part models it is decomposed exactly into extensive (participation / non-structural-zero) and intensive (count) channels that sum to the total.

Usage

## S3 method for class 'gec'
first_difference(object, variable, from, to, level = 0.95, ...)

## S3 method for class 'gec_fe'
first_difference(object, variable, from, to, level = 0.95, ...)

## S3 method for class 'cpb_fe'
first_difference(object, variable, from, to, level = 0.95, ...)

## S3 method for class 'hurdle_gec'
first_difference(
  object,
  variable,
  from,
  to,
  level = 0.95,
  stage = c("both", "participation", "intensity", "zero", "count"),
  ...
)

## S3 method for class 'zi_gec'
first_difference(
  object,
  variable,
  from,
  to,
  level = 0.95,
  stage = c("both", "participation", "intensity", "zero", "count"),
  ...
)

## S3 method for class 'count_reg'
first_difference(object, variable, from, to, level = 0.95, ...)

## S3 method for class 'hurdle_count'
first_difference(
  object,
  variable,
  from,
  to,
  level = 0.95,
  stage = c("both", "participation", "intensity", "zero", "count"),
  ...
)

## S3 method for class 'zi_count'
first_difference(
  object,
  variable,
  from,
  to,
  level = 0.95,
  stage = c("both", "participation", "intensity", "zero", "count"),
  ...
)

Arguments

Value

A "ud_fd" data frame – the package-wide first-difference contract shared by every first_difference() method: columns component, from, to, diff, lower, upper, method. Single-equation fits return one row (component = "mean"); two-part fits return one row per margin level (binary-stage probability, count-stage mean, marginal mean). method records the uncertainty source per row ("delta", a bootstrap label, or "none" with NA bounds).

Extensive/intensive first-difference decomposition for a hurdle- or ZI-CPB

Description

The effect of moving variable from one value to another, decomposed into the change in the participation (extensive) margin, the intensity (intensive) margin, and the marginal expected count, holding other covariates at their means. This is the quantity the single-equation defaults cannot separate.

Usage

## S3 method for class 'hurdle_cpb'
first_difference(
  object,
  variable,
  from,
  to,
  B = 0,
  level = 0.95,
  data = NULL,
  stage = c("both", "participation", "intensity", "zero", "count"),
  ...
)

## S3 method for class 'zi_cpb'
first_difference(
  object,
  variable,
  from,
  to,
  B = 0,
  level = 0.95,
  data = NULL,
  stage = c("both", "participation", "intensity", "zero", "count"),
  ...
)

Arguments

Value

A "ud_fd" data frame (the package-wide first-difference contract): columns component, from, to, diff, lower, upper, method, one row per margin level.

Fitted means for the two-part and fixed-effects classes

Description

Completes stats::fitted() across the family, returning the same quantity as predict(object, type = "response") on the estimation data: the mean for the single-equation fits and the marginal expected count (zeros included) for the two-part fits. gec_fe inherits fitted.gec.

Usage

## S3 method for class 'cpb_fe'
fitted(object, ...)

## S3 method for class 'hurdle_cpb'
fitted(object, ...)

## S3 method for class 'hurdle_gec'
fitted(object, ...)

## S3 method for class 'hurdle_count'
fitted(object, ...)

## S3 method for class 'zi_cpb'
fitted(object, ...)

## S3 method for class 'zi_gec'
fitted(object, ...)

## S3 method for class 'zi_count'
fitted(object, ...)

Arguments

Value

A numeric vector, one element per estimation observation.

Generalized event count (Katz-family) regression

Description

Fits King's generalized event count model, a Katz-family count regression whose single dispersion parameter delta (the variance-to-mean ratio) is estimated freely and spans underdispersion (delta < 1, a finite-support member; the continuous parameter binomial is this cell), equidispersion (delta = 1, Poisson), and overdispersion (delta > 1, the negative binomial). Unlike cpb(), which fixes the direction of dispersion to under, gec() lets the data choose. The Katz recursion delivers exact first and second moments—the Winkelmann–Signorino–King correction realized directly—and the likelihood is evaluated in C++.

Usage

gec(
  formula,
  data,
  truncated = FALSE,
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  offset = NULL,
  max.support = 500,
  maxit = 20000,
  reltol = 1e-08
)

Arguments

Value

An object of class "gec" with coefficients, delta (the estimated dispersion), loglik, bootstrap standard errors, and bookkeeping.

See Also

cpb(), compare_dispersion()

Examples

set.seed(1); x <- rnorm(400)
y <- rpois(400, exp(1 + 0.5 * x))
gec(y ~ x, data = data.frame(y = y, x = x), se = "none")

Generalized event count (Katz) distribution functions

Description

Density, distribution, quantile, and random generation for the King (1989) generalized event count model with rate lambda (the mean) and dispersion delta (the variance-to-mean ratio: ⁠< 1⁠ underdispersed, ⁠= 1⁠ Poisson, ⁠> 1⁠ overdispersed). Consistent with the estimator gec().

Usage

dgec(x, lambda, delta, max.support = 500, log = FALSE)

pgec(q, lambda, delta, max.support = 500, lower.tail = TRUE, log.p = FALSE)

qgec(p, lambda, delta, max.support = 500, lower.tail = TRUE, log.p = FALSE)

rgec(n, lambda, delta, max.support = 500)

Arguments

Value

dgec a density, pgec a CDF, qgec a quantile, rgec a numeric vector of count draws.

See Also

gec()

Examples

dgec(0:5, lambda = 3, delta = 0.7)
var(rgec(2000, lambda = 3, delta = 0.7)) / mean(rgec(2000, lambda = 3, delta = 0.7))

GEC (Katz-family) regression with high-dimensional unit fixed effects

Description

Fits gec() with a full set of unit fixed effects, concentrating (profiling) out the unit intercepts by a one-dimensional inner maximization per unit, so the outer optimizer handles only the covariate coefficients and the free dispersion parameter delta. This is the cpb_fe() concentration generalized to the whole Katz family: the panel can be under-, equi-, or overdispersed and the direction is estimated, not presumed.

Usage

gec_fe(
  formula,
  data,
  fe,
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  offset = NULL,
  max.support = 500L,
  inner_it = 30L,
  maxit = 3000L,
  reltol = 1e-07,
  bias_correct = c("none", "jackknife")
)

Arguments

Details

Limitation: unlike cpb_fe(), gec_fe() has no truncated argument – a concentrated zero-truncated GEC is not currently implemented. For a zero-truncated GEC with fixed effects, enter the unit factor as dummies in gec()'s formula (feasible for moderate unit counts).

Short panels: delta and the fixed effects carry the incidental-parameters bias of nonlinear fixed-effects estimation, of order 1/T; treat delta cautiously below about T = 30. The covariate coefficients are not materially affected. bias_correct = "jackknife" removes the leading 1/T term by the split-panel jackknife exactly as in cpb_fe(): refit on each unit's temporal halves and report 2 * full - mean(halves), with the unit effects and fitted values re-concentrated at the corrected parameters and logLik/AIC kept at the maximum-likelihood fit (uncorrected estimates in ⁠$uncorrected⁠). The same time-homogeneity validity gate applies: on panels where the two halves do not estimate a common parameter the correction is REFUSED with a warning and the maximum-likelihood fit is returned (see cpb_fe(), Details).

Value

An object of class c("gec_fe", "gec").

See Also

gec(), cpb_fe()

Examples

set.seed(5)
d <- do.call(rbind, lapply(1:30, function(i) {
  x <- rnorm(10); N <- pmax(round(exp(1 + rnorm(1, 0, 0.4) + 0.3 * x) / 0.5), 1)
  data.frame(unit = i, x = x, y = rbinom(10, N, 0.5))
}))
gec_fe(y ~ x, data = d, fe = "unit")   # delta ~ 0.5: underdispersed

Glance at a CPB fit (broom method)

Description

Glance at a CPB fit (broom method)

Usage

## S3 method for class 'cpb'
glance(x, ...)

## S3 method for class 'cpb_fe'
glance(x, ...)

## S3 method for class 'hurdle_cpb'
glance(x, ...)

## S3 method for class 'zi_cpb'
glance(x, ...)

## S3 method for class 'gec'
glance(x, ...)

## S3 method for class 'hurdle_gec'
glance(x, ...)

## S3 method for class 'zi_gec'
glance(x, ...)

## S3 method for class 'count_reg'
glance(x, ...)

## S3 method for class 'zi_count'
glance(x, ...)

## S3 method for class 'hurdle_count'
glance(x, ...)

Arguments

Value

A one-row data frame of fit statistics.

Hurdle Poisson / negative-binomial regression (matched CPB baseline)

Description

Fits a participation logit joined to a zero-truncated Poisson or negative-binomial intensity, the count analogue of hurdle_cpb(). Returned as a "hurdle_count" object that compare_models() and score() accept, so a hurdle-CPB and a hurdle-NB can be compared on one footing.

Usage

hurdle_count(
  formula,
  data,
  family = c("poisson", "negbin", "compois"),
  participation = NULL,
  fe = NULL,
  part_fe = NULL,
  link = c("logit", "probit", "cloglog"),
  offset = NULL,
  se = c("analytic", "robust", "cluster", "none"),
  cluster = NULL,
  ...
)

Arguments

Value

An object of class "hurdle_count".

See Also

hurdle_cpb(), zi_count(), compare_models()

Examples

set.seed(1); n <- 300; x <- rnorm(n); z <- rnorm(n)
y <- ifelse(rbinom(n, 1, plogis(0.4 + 0.8 * z)) == 1, rpois(n, exp(1 + 0.3 * x)) + 1L, 0L)
hurdle_count(y ~ x, data.frame(y = y, x = x, z = z), family = "poisson",
             participation = ~ z, se = "none")

Hurdle continuous parameter binomial regression

Description

Fits a participation (hurdle) model joined to a zero-truncated CPB for the positive counts. Because the hurdle log-likelihood factorizes, the two parts are fit separately: a logistic regression of participation on all units, and a zero-truncated CPB (cpb()) on the positive counts. The result is the natural model for a bounded, underdispersed participation process — most units at zero, participants carrying a tight count.

Usage

hurdle_cpb(
  formula,
  data,
  participation = NULL,
  fe = NULL,
  part_fe = NULL,
  link = c("logit", "probit", "cloglog"),
  offset = NULL,
  cluster = NULL,
  se = c("none", "bootstrap"),
  B = 500,
  ...
)

Arguments

Value

An object of class "hurdle_cpb" with elements participation (a glm), intensity (a cpb or cpb_fe), and bookkeeping.

Examples

set.seed(1)
n <- 800; x <- rnorm(n); z <- rnorm(n)
y <- rhurdle_cpb(n, lambda = exp(1.3 + 0.5 * x), alpha = 0.5,
                 p = plogis(-0.2 + 0.8 * z))
fit <- hurdle_cpb(y ~ x, data = data.frame(y = y, x = x, z = z),
                  participation = ~ z)
fit

Hurdle GEC (Katz-family) regression

Description

Joins a participation (hurdle) logit to a zero-truncated GEC intensity on the positive counts. Unlike hurdle_cpb(), whose intensity is fixed to the underdispersed CPB, the GEC intensity estimates its own dispersion direction.

Usage

hurdle_gec(
  formula,
  data,
  participation = NULL,
  part_fe = NULL,
  link = c("logit", "probit", "cloglog"),
  offset = NULL,
  cluster = NULL,
  se = c("none", "bootstrap"),
  B = 500,
  max.support = 500
)

Arguments

Value

An object of class "hurdle_gec".

Limitation

Unlike hurdle_cpb(), hurdle_gec() has no fe argument for a concentrated fixed-effects intensity (no concentrated zero-truncated GEC is implemented). For within-unit intensities, enter the unit factor as dummies in formula (feasible for moderate unit counts), or use hurdle_cpb().

See Also

hurdle_cpb(), gec(), zi_gec()

Examples

set.seed(3); n <- 300; x <- rnorm(n); z <- rnorm(n)
y <- ifelse(rbinom(n, 1, plogis(0.3 + 0.8 * z)) == 1, rpois(n, exp(1 + 0.3 * x)) + 1L, 0L)
hurdle_gec(y ~ x, data.frame(y = y, x = x, z = z), participation = ~ z)

Implied ceiling with a profile-likelihood interval

Description

Returns the observation- (or profile-) specific ceiling lambda/(1-alpha), with an interval propagating the profile-likelihood uncertainty in alpha at the fitted mean. (Coefficient uncertainty in lambda is not propagated here; use first_difference() with quantity = "ceiling" for a fully bootstrapped contrast.)

Usage

implied_ceiling(object, ...)

## S3 method for class 'cpb'
implied_ceiling(object, newdata = NULL, level = 0.95, ...)

Arguments

Value

A data frame with lambda, ceiling, and lower/upper bounds.

Examples

set.seed(6); x <- rnorm(300)
N <- pmax(round(exp(1.5 + 0.4 * x) / 0.5), 1); y <- rbinom(300, N, 0.5)
fit <- cpb(y ~ x, data.frame(y = y, x = x)[y > 0, ], se = "none")
implied_ceiling(fit, newdata = data.frame(x = c(-1, 0, 1)))

Incidence rate ratios for a CPB fit

Description

Incidence rate ratios for a CPB fit

Usage

irr(object, ...)

## S3 method for class 'cpb'
irr(object, level = 0.95, ...)

Arguments

Value

A "ud_irr" data frame – the package-wide rate-ratio contract (columns term, equation, ratio, estimate, lower, upper, method) shared by every irr() method; here with bootstrap percentile intervals from the stored draws.

Examples

set.seed(8); x <- rnorm(300)
N <- pmax(round(exp(1.5 + 0.4 * x) / 0.5), 1); y <- rbinom(300, N, 0.5)
fit <- cpb(y ~ x, data.frame(y = y, x = x)[y > 0, ], se = "bootstrap", B = 100)
irr(fit)

Rate and odds ratios for the matched count families

Description

Rate ratios (exp of the count/intensity coefficients, column IRR) and, for the two-part models, odds ratios of the binary stage (column OR), by equation. Element names match the CPB family (intensity for the count component; participation for the hurdle stage, zero for the inflation stage).

Usage

## S3 method for class 'count_reg'
irr(object, level = 0.95, ...)

## S3 method for class 'gec'
irr(object, level = 0.95, ...)

## S3 method for class 'gec_fe'
irr(object, level = 0.95, ...)

## S3 method for class 'cpb_fe'
irr(object, level = 0.95, ...)

## S3 method for class 'hurdle_gec'
irr(object, level = 0.95, ...)

## S3 method for class 'zi_gec'
irr(object, level = 0.95, ...)

## S3 method for class 'hurdle_count'
irr(object, level = 0.95, ...)

## S3 method for class 'zi_count'
irr(object, level = 0.95, ...)

Arguments

Value

A "ud_irr" data frame – the package-wide ratio contract shared by every irr() method: columns term, equation ("count" or "binary"), ratio ("IRR" or "OR"), estimate, lower, upper, and method (the interval source; "none" with NA bounds when no covariance is available). Two-part models stack both equations.

Mundlak (correlated random effects) device

Description

Augments a data frame with the unit-level means of the time-varying numeric covariates in formula, and returns the augmented formula and data. Fitting any of the package's estimators on the result implements the Mundlak / correlated-random-effects specification: the coefficients on the original covariates recover the within-unit (fixed-effects-consistent) effects, while the coefficients on the unit means capture (and test) the correlation between the covariates and the unit effect. It is a lighter, between-variation-preserving alternative to full fixed effects, and it is model-agnostic — the same device feeds cpb(), gec(), or count_reg().

Usage

mundlak(formula, data, unit, suffix = "_mean")

Arguments

Value

A list with formula (the augmented formula), data (the augmented data frame), and added (the names of the unit-mean columns).

See Also

count_reg(), cpb()

Examples

set.seed(1)
d <- data.frame(y = rpois(200, 3), x = rnorm(200), unit = factor(rep(1:20, each = 10)))
m <- mundlak(y ~ x, d, unit = "unit")
fit <- count_reg(m$formula, m$data, family = "poisson")

UN peacekeeping contributions, state-years

Description

A state-year panel of contributions to United Nations peacekeeping operations, the archetype of a bounded, underdispersed participation count: most state-years contribute to no operation, and the states that contribute carry a tight count held in a narrow band by a turning-over roster. About 42\ hurdle-CPB.

Usage

data(peacekeeping)

Format

A data frame with 4,448 rows and 9 variables:

iso3

ISO3 country code (the panel unit).

year

Calendar year, 1990–2024.

contributions

Number of distinct UN peacekeeping operations the state contributes personnel to in that year (the count outcome).

democracy

Electoral-democracy index.

lgdppc

Log GDP per capita.

lpop

Log population.

milper

Military personnel, thousands.

majorpower

Major-power indicator (0/1).

region

World region.

Source

Contribution counts aggregated from the International Peace Institute Providing for Peacekeeping database (https://www.providingforpeacekeeping.org); covariates from V-Dem and the Correlates of War National Material Capabilities data. Counts and public covariates only; prepared for illustration.

Examples

data(peacekeeping)
table(peacekeeping$contributions == 0)

fit <- hurdle_cpb(contributions ~ lgdppc + milper, data = peacekeeping,
                  participation = ~ democracy + majorpower, fe = "iso3")
fit

Non-randomized PIT histogram for a fitted count model

Description

Draws the non-randomized probability integral transform histogram of Czado, Gneiting, and Held (2009). A well-calibrated model yields a flat histogram at height one (the reference line); a U shape indicates under-dispersion in the predictive distribution and a hump indicates over-dispersion.

Usage

pit_hist(
  fit,
  bins = 10,
  main = "PIT histogram",
  xlab = "PIT",
  ylab = "Relative frequency",
  ...
)

Arguments

Value

Invisibly, the vector of bin heights (normalized so that a calibrated model gives heights near one).

Examples

set.seed(1)
y <- rcpb(400, lambda = 3, alpha = 0.5)
fit <- cpb(y ~ 1, data = data.frame(y = y), truncated = FALSE, se = "none")
pit_hist(fit)

Predictions from a matched Poisson/NB/COM-Poisson fit

Description

Predictions from a matched Poisson/NB/COM-Poisson fit

Usage

## S3 method for class 'count_reg'
predict(
  object,
  newdata = NULL,
  type = c("response", "link", "prob"),
  at = NULL,
  offset = NULL,
  ...
)

Arguments

Value

A numeric vector.

Predictions from a CPB fit

Description

Predictions from a CPB fit

Usage

## S3 method for class 'cpb'
predict(
  object,
  newdata = NULL,
  type = c("response", "link", "ceiling", "prob"),
  at = NULL,
  offset = NULL,
  ...
)

Arguments

Value

A numeric vector.

Examples

set.seed(1); x <- rnorm(300)
N <- pmax(round(exp(1.6 + 0.5 * x) / 0.5), 1); y <- rbinom(300, N, 0.5)
fit <- cpb(y ~ x, data = data.frame(y = y, x = x)[y > 0, ], se = "none")
predict(fit, newdata = data.frame(x = 0), type = "ceiling")
predict(fit, newdata = data.frame(x = c(-1, 1)), type = "prob", at = 5)

Predictions from a fixed-effects CPB fit

Description

Predictions from a fixed-effects CPB fit

Usage

## S3 method for class 'cpb_fe'
predict(object, newdata = NULL, type = c("response", "link", "ceiling"), ...)

Arguments

Value

A numeric vector.

Predictions from a generalized event count fit

Description

Predictions from a generalized event count fit

Usage

## S3 method for class 'gec'
predict(
  object,
  newdata = NULL,
  type = c("response", "link", "prob"),
  at = NULL,
  offset = NULL,
  ...
)

Arguments

Value

A numeric vector.

Predict from a hurdle Poisson/NB fit

Description

Predict from a hurdle Poisson/NB fit

Usage

## S3 method for class 'hurdle_count'
predict(
  object,
  newdata = NULL,
  type = c("response", "participation", "intensity"),
  offset = NULL,
  ...
)

Arguments

Value

A numeric vector.

Predict from a hurdle-CPB fit

Description

Predict from a hurdle-CPB fit

Usage

## S3 method for class 'hurdle_cpb'
predict(
  object,
  newdata = NULL,
  type = c("response", "participation", "intensity"),
  ...
)

Arguments

Value

A numeric vector.

Predict from a hurdle-GEC fit

Description

Predict from a hurdle-GEC fit

Usage

## S3 method for class 'hurdle_gec'
predict(
  object,
  newdata = NULL,
  type = c("response", "participation", "intensity"),
  ...
)

Arguments

Value

A numeric vector.

Predict from a zero-inflated Poisson/NB fit

Description

Predict from a zero-inflated Poisson/NB fit

Usage

## S3 method for class 'zi_count'
predict(
  object,
  newdata = NULL,
  type = c("response", "count", "zero", "intensity"),
  offset = NULL,
  ...
)

Arguments

Value

A numeric vector.

Predict from a zi_cpb fit

Description

Predict from a zi_cpb fit

Usage

## S3 method for class 'zi_cpb'
predict(object, newdata = NULL, type = c("response", "zero", "intensity"), ...)

Arguments

Value

A numeric vector.

Predict from a zi_gec fit

Description

Predict from a zi_gec fit

Usage

## S3 method for class 'zi_gec'
predict(object, newdata = NULL, type = c("response", "zero", "intensity"), ...)

Arguments

Value

A numeric vector.

Simulate from the continuous parameter binomial

Description

Draws counts from the continuous parameter binomial (CPB) or its zero-truncated variant.

Usage

rcpb(n, lambda, alpha, truncated = FALSE)

Arguments

Value

An integer vector of counts.

Examples

set.seed(1)
table(rcpb(1000, lambda = 3, alpha = 0.5))

Simulate from the hurdle continuous parameter binomial

Description

Simulate from the hurdle continuous parameter binomial

Usage

rhurdle_cpb(n, lambda, alpha, p)

Arguments

Value

An integer vector: 0 for nonparticipants, a zero-truncated CPB draw otherwise.

Examples

set.seed(1)
y <- rhurdle_cpb(1000, lambda = 3, alpha = 0.5, p = 0.6)
mean(y == 0)

Hanging rootogram for a fitted count model

Description

Draws a Tukey hanging rootogram: bars for the observed frequencies hang from the curve of expected frequencies, both on the square-root scale. Bars that hang below the zero line mark counts the model under-predicts; bars that stop short mark counts it over-predicts.

Usage

rootogram(
  fit,
  kmax = NULL,
  main = "Rootogram",
  xlab = "Count",
  ylab = "sqrt(frequency)",
  ...
)

Arguments

Value

Invisibly, a data frame of count, observed, and expected frequencies.

Examples

set.seed(1)
y <- rcpb(400, lambda = 3, alpha = 0.5)
fit <- cpb(y ~ 1, data = data.frame(y = y), truncated = FALSE, se = "none")
rootogram(fit)

Simulate from the zero-inflated continuous parameter binomial

Description

Simulate from the zero-inflated continuous parameter binomial

Usage

rzicpb(n, lambda, alpha, pi)

Arguments

Value

An integer vector: a structural zero with probability pi, otherwise an (untruncated) CPB draw, which may itself be zero.

Examples

set.seed(1)
y <- rzicpb(1000, lambda = 3, alpha = 0.5, pi = 0.3)
mean(y == 0)

Proper scoring rules for a fitted count model

Description

Computes the mean logarithmic score and the ranked probability score (RPS) of a fitted model's predicted distribution against observed counts. Lower is better for both, and both are proper, so they are the natural way to compare the calibration of competing count models.

Usage

score(fit, newdata = NULL, kmax = NULL)

Arguments

Details

By default the scores are computed in sample (against the data the model was fit to). Supply newdata to score a fitted model on held-out observations, or use cv_score() for a cross-validated score; an in-sample log score equals -logLik/n and does not penalize model complexity, so for comparing models of different size the held-out or cross-validated score is the honest criterion.

Value

A named numeric vector c(logscore, rps).

See Also

cv_score()

Examples

set.seed(1)
y <- rcpb(400, lambda = 3, alpha = 0.5)
fit <- cpb(y ~ 1, data = data.frame(y = y), truncated = FALSE, se = "none")
score(fit)

Simulate responses from a fitted underdisp model

Description

Draws nsim replicate response vectors from the fitted model, at the estimated parameters and the estimation data, in the format of stats::simulate(). Available for every model class in the package: cpb, cpb_fe, gec, gec_fe (through inheritance), count_reg, hurdle_count, zi_count, hurdle_cpb, hurdle_gec, zi_cpb, and zi_gec. Two-part models first draw the binary stage (participation or structural zero), then the count stage from its own distribution, so the replicates carry the model's full zero structure.

Usage

## S3 method for class 'cpb'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'cpb_fe'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'gec'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'count_reg'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'hurdle_cpb'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'hurdle_gec'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'hurdle_count'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'zi_cpb'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'zi_gec'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'zi_count'
simulate(object, nsim = 1, seed = NULL, ...)

Arguments

Details

The main consumer is simulated-residual diagnostics: DHARMa::createDHARMa(simulatedResponse = as.matrix(simulate(fit, 250)), observedResponse = y, fittedPredictedResponse = fitted(fit), integerResponse = TRUE) works for any fit in the family.

Value

A data frame with nsim integer columns, one row per observation, with a "seed" attribute.

Examples

set.seed(1); x <- rnorm(200)
N <- pmax(round(exp(1.4 + 0.4 * x) / 0.5), 1); y <- rbinom(200, N, 0.5)
fit <- cpb(y ~ x, data = data.frame(y = y, x = x)[y > 0, ], se = "none")
sims <- simulate(fit, nsim = 5)
colMeans(sims)

Summarize a CPB fit

Description

Summarize a CPB fit

Usage

## S3 method for class 'cpb'
summary(object, ...)

Arguments

Value

An object of class "summary.cpb" with the coefficient table, the dispersion parameter and its profile-likelihood interval, the implied ceiling, fit statistics, and the likelihood-ratio test against a (zero-truncated) Poisson.

Tidy a CPB fit (broom method)

Description

Tidy a CPB fit (broom method)

Usage

## S3 method for class 'cpb'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'cpb_fe'
tidy(x, ...)

## S3 method for class 'hurdle_cpb'
tidy(x, ...)

## S3 method for class 'zi_cpb'
tidy(x, ...)

## S3 method for class 'gec'
tidy(x, ...)

## S3 method for class 'hurdle_gec'
tidy(x, ...)

## S3 method for class 'zi_gec'
tidy(x, ...)

## S3 method for class 'count_reg'
tidy(x, ...)

## S3 method for class 'zi_count'
tidy(x, ...)

## S3 method for class 'hurdle_count'
tidy(x, ...)

Arguments

Value

A data frame with one row per coefficient.

Screen a count outcome for underdispersion

Description

Applies the diagnostic sequence developed in Bagozzi (2026): a marginal verdict from the conditional Pearson statistic and a through-origin score regression (which gives the direction of dispersion), a negative-binomial-versus-Poisson test for the overdispersion call, and—most importantly—an at-risk verdict on the positive counts benchmarked against a zero-truncated Poisson. The last step is what separates genuine underdispersion from the artifact of conditioning on Y>0. When the data are not zero-dominated it also fits the CPB and reports a ceiling-exceedance diagnostic, and it fits a generalized-Poisson soft-tail comparator.

Usage

ud_screen(
  formula,
  data,
  run_cpb = TRUE,
  cpb_max_n = 3000,
  run_gp = TRUE,
  run_comp = TRUE,
  comp_max_par = 30,
  comp_max_n = 5000,
  ztp_threshold = c("calibrated", "bootstrap"),
  ztp_boot_B = 199L,
  digits = 3
)

Arguments

Value

An object of class "ud_screen" with verdict_marginal, verdict_atrisk, the conditional and at-risk (ZTP-benchmarked) Pearson statistics, the NB-vs-Poisson LR test, a log-likelihood comparison, (when fit) the CPB alpha and ceiling-exceedance share, and the over-conditioning guard state: atrisk_skipped (TRUE when the mean model nearly saturates the positive counts, so the at-risk statistic is not computed and the printout says why), sat_ratio (the fitted parameter share of the positives), and overconditioned (TRUE when that share reaches 0.10, the region where the calibrated threshold is anti-conservative; the printout then flags the verdict as diagnostic rather than probative and recommends ztp_threshold = "bootstrap").

References

King, G. (1989). Variance specification in event count models. AJPS 33(3):762-784.

Bagozzi, B. E. (2026). Revisiting underdispersion in political science. Companion manuscript.

See Also

cpb(), compare_dispersion()

Examples

set.seed(1); x <- rnorm(250)
N <- pmax(round(exp(1.4 + 0.4 * x) / 0.6), 1); y <- rbinom(250, N, 0.6)
ud_screen(y ~ x, data = data.frame(y = y, x = x))

Bootstrap covariance for a fixed-effects CPB fit

Description

The covariance of the covariate coefficients from the stored pairs/cluster bootstrap (se = "bootstrap" at fit time); NULL when no bootstrap was run.

Usage

## S3 method for class 'cpb_fe'
vcov(object, ...)

Arguments

Value

A covariance matrix, or NULL.

Zero-inflated Poisson / negative-binomial regression (matched CPB baseline)

Description

Fits a structural-zero mixture with a Poisson or negative-binomial count component, the count analogue of zi_cpb(). The count component uses a log link (on the mean for Poisson/NB, on the rate \lambda for COM-Poisson) and the inflation probability a link set by link. Returned as a "zi_count" object that compare_models() and score() accept.

Usage

zi_count(
  formula,
  data,
  family = c("poisson", "negbin", "compois"),
  zero = NULL,
  fe = NULL,
  zero_fe = NULL,
  link = c("logit", "probit", "cloglog"),
  offset = NULL,
  se = c("analytic", "robust", "cluster", "none"),
  cluster = NULL,
  ...
)

Arguments

Value

An object of class "zi_count".

See Also

zi_cpb(), hurdle_count(), compare_models()

Examples

set.seed(2); n <- 300; x <- rnorm(n); z <- rnorm(n)
y <- ifelse(rbinom(n, 1, plogis(-0.5 + 0.8 * z)) == 1, 0L, rpois(n, exp(1 + 0.3 * x)))
zi_count(y ~ x, data.frame(y = y, x = x, z = z), family = "poisson",
         zero = ~ z, se = "none")

Zero-inflated continuous parameter binomial regression

Description

Fits the zero-inflated CPB. The model mixes a structural-zero process (a logistic model for the probability pi that a unit is a structural zero) with an untruncated CPB for the count, so a zero can arise either structurally or as a sampling zero from the CPB. Unlike hurdle_cpb(), the likelihood does not factorize, so the parameters are estimated by direct joint maximum likelihood, seeded from separate fits (a zero-truncated CPB on the positives and a logit of the zero indicator) and cross-checked against an expectation-maximization climb of the same observed-data likelihood, keeping the better optimum. An EM-only path (method = "em") is retained for comparison.

Usage

zi_cpb(
  formula,
  data,
  zero = NULL,
  fe = NULL,
  zero_fe = NULL,
  method = c("ml", "em"),
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  max.support = 500,
  maxit = 200,
  tol = 1e-06,
  offset = NULL
)

Arguments

Details

The inflation equation uses a logit link. For a probit or cloglog inflation link, use zi_count() (Poisson/NB/COM-Poisson count), whose link argument covers the binary stage.

Value

An object of class "zi_cpb".

Examples

set.seed(1)
n <- 800; x <- rnorm(n); z <- rnorm(n)
y <- rzicpb(n, lambda = exp(1.3 + 0.5 * x), alpha = 0.5, pi = plogis(-0.5 + 0.8 * z))
zi_cpb(y ~ x, data = data.frame(y = y, x = x, z = z), zero = ~ z)

Zero-inflated GEC (Katz-family) regression

Description

A structural-zero mixture with a GEC count component whose dispersion delta is estimated freely. Estimated by robust direct maximum likelihood, seeded from separate fits (a GEC on the positive counts for the count parameters and a logit of the zero indicator for the inflation), exactly as in zi_cpb(); this avoids the degenerate pi -> 0 basin that coordinate-wise EM falls into.

Usage

zi_gec(
  formula,
  data,
  zero = NULL,
  zero_fe = NULL,
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  max.support = 500,
  maxit = 200,
  tol = 1e-06,
  offset = NULL
)

Arguments

Value

An object of class "zi_gec".

See Also

zi_cpb(), gec(), hurdle_gec()

Examples

set.seed(4); n <- 300; x <- rnorm(n); z <- rnorm(n)
y <- ifelse(rbinom(n, 1, plogis(-0.5 + 0.8 * z)) == 1, 0L, rpois(n, exp(1 + 0.3 * x)))
zi_gec(y ~ x, data.frame(y = y, x = x, z = z), zero = ~ z)

Boundary-corrected test for zero-inflation (CPB vs ZI-CPB)

Description

A likelihood-ratio test of whether a count needs a structural-zero component. The (untruncated) CPB is nested in the zero-inflated CPB at a structural-zero probability of zero. Because that value lies on the boundary of the parameter space, the LR statistic follows a \tfrac12\chi^2_0 + \tfrac12\chi^2_1 mixture (Self and Liang 1987), which halves the naive \chi^2_1 p-value. The test uses an intercept-only structural-zero probability, so it is a clean single-parameter boundary test; a covariate-dependent structural-zero model is better compared with information criteria and proper scores via compare_dispersion().

Usage

zi_test(object, object2 = NULL, data = NULL, ...)

Arguments

Value

An object of class "zi_test" with the two log-likelihoods, the LR statistic, and the boundary-corrected p-value.

Examples

set.seed(1); x <- rnorm(500)
y <- rzicpb(500, lambda = exp(1.2 + 0.4 * x), alpha = 0.5, pi = 0.3)
zi_test(y ~ x, data = data.frame(y = y, x = x))