Changelog • underdisp

underdisp 0.2.1

Tests

A test that dispersion_test() gives the same statistic for a cpb_fe() fit whether or not a covariate is centred allowed a relative difference of 1e-6 and failed on CRAN’s Linux ARM64 check, where the two statistics differed by 1.6e-6 (4e-7 on x86_64). The gap is where the optimizer of the fixed-effects fit stops; the test now allows 1e-4. It guards against the error fixed in 0.2.0, which on the test’s data gave the uncentred covariate a statistic of 2e10. The package code is unchanged.

Documentation

?ud_screen described the constant of the calibrated at-risk threshold, 1 - 2.27/sqrt(n+), as an estimated standard deviation. It is a scale matched to the 5% quantile of the simulated null on one calibration design, where the null is right-skewed and its standard deviation somewhat larger. The threshold is unchanged.

underdisp 0.2.0

CRAN release: 2026-11-04

The continuous parameter binomial (CPB) now follows King’s (1989) definition by default; support = "floor" reproduces the CPB fits of 0.1.x. The release adds two families, average marginal effects, prediction intervals and separation checks, and fixes errors in fixed-effects, two-part and new-data computations, some of which change earlier results.

Breaking changes

The CPB follows King’s definition

Every CPB function defaults to King’s CPB as Winkelmann, Signorino and King
1. normalize it, on 0, ..., ceiling(lambda / (1 - alpha)): the generalized event count of gec() at delta = alpha < 1. Earlier versions used 0, ..., floor(lambda / (1 - alpha)), a different distribution; the two agree when the ceiling is an integer. King’s CPB has a mean at or above lambda and a continuous log-likelihood that tends to the Poisson’s as alpha approaches 1, so max.support no longer bounds alpha; it caps only the width of the pmf the likelihood evaluates.
support = "floor" gives the earlier definition in cpb(), cpb_fe(), hurdle_cpb(), zi_cpb(), dcpb(), pcpb(), qcpb(), rcpb(), rhurdle_cpb(), rzicpb(), compare_dispersion() and ud_screen(), and reproduces the estimates, log-likelihoods and degrees of freedom of 0.1.x (zero-inflated fits can reach higher maxima from the new starting values). Fixes that do not depend on the support, such as those to intervals and to the reference unit, apply to it too. Objects saved by earlier versions are read as floor-support fits.
To reproduce 0.1.x, the floor support keeps the earlier behavior of: the cpb_fe() sentinel for units with no positive count (also in dispersion_test()); the profile trace for alpha (down to 0.005); the zi_cpb(fe =) plug-in fit and filled-in unit effects; zi_cpb()’s count mean (the rate) and its EM convergence flag; in hurdle_cpb() and zi_cpb(), the degrees of freedom and new-data predictions of units at a limit and the held-out scoring of unseen units; uncorrected replicates of jackknife-corrected fits; and max.support = 500 in zi_cpb() and compare_dispersion().
On King’s support a fit with alpha above 0.999 warns that the data show no underdispersion (dispersion_test() and summary() then report a statistic of 0), and a fit with alpha below 0.001 notes that the interval and tests for alpha are at the two-point boundary.
On King’s support the profile interval for alpha is traced down to 1e-6 (a lower side that never reaches its cut ends at 0); predict() and fitted() of zi_cpb() use the exact mean of the count component; and the distribution functions, predict(type = "prob"), simulate() and calibrate_alpha() return NA with a warning past the 1e5-term cap.
The documentation reports simulation results on King’s support: the first-order profile interval for alpha covers 0.940 to 0.956 (0.90 to 0.94 on the floor support, where calibrate_alpha() is the remedy); the bootstrap’s 95% slope interval covers 0.91 to 0.95 at n = 300 (0.95 on the floor support); cpb_fe()’s bias in alpha at T = 6 is -0.122, and +0.008 after the jackknife (floor support: -0.153).

Other changes to defaults and results

max.support = NULL, meaning max(500, 10 * max(y)), is the default of gec(), hurdle_gec(), zi_gec(), zi_cpb() and compare_dispersion(), as it was of cpb(), cpb_fe() and gec_fe() (with the old 500, gec() failed on counts far above it).
bias_correct = "jackknife" requires time =, the column that orders each unit’s rows; rows not in time order were split as they came, which defeated the validity gate. With se = "bootstrap" every replicate is corrected too, so the standard errors include the correction’s variance.
Every first_difference() method takes its interval by one rule: the percentile interval over the stored bootstrap draws, jointly over the coefficients, the dispersion parameter and a fixed-effects fit’s reference unit (draws recentered on the estimate), or else the delta method over the full covariance. The intervals of gec(), gec_fe(), cpb_fe() and count_reg() held the dispersion parameter fixed, and some two-part rows had none.
The hurdle and zero-inflated fits do not count in their degrees of freedom the parameters of a unit with no positive count, which sit at their limits. AIC and BIC change on panels with such units: the hurdles counted their participation dummies while zi_cpb(fe =) did not, favoring the mixture by two AIC points per unit, and zi_count() counted two parameters per unit. See ?compare_models.
zi_cpb(fe =) leaves units with no positive count out of the estimation and says how many (they used to get a filled-in effect whose zeros shaped the inflation estimates). With method = "ml" it polishes its plug-in fit with every unit’s effect free, so it is a joint maximum; the plug-in could lie tens of log-likelihood points lower and depended on the covariates’ centering.
zi_cpb(method = "em") reports the observed-data log-likelihood at its estimate (a floor of 1e-12 on each probability overstated it) and, on King’s support, whether it converged.
predict(newdata) and implied_ceiling(newdata) on cpb_fe(), gec_fe() and zi_cpb(fe =) fits enter a unit the fit estimated at its own effect, not at the reference unit. These fits return fitted(), residuals(), predict() and simulate() in the data’s row order, not sorted by unit.
A fit with an offset carries it to new rows: predict(), score(), cv_score() and ceiling_exceedance() read it from newdata or take offset =, and refuse rows without one. They used exposure 1.
implied_ceiling()’s lower and upper are renamed ceiling_at_alpha_lower and ceiling_at_alpha_upper (they are not a confidence interval), and the new max_count is the largest count with positive probability.
compare_dispersion() and ud_screen() drop the in-sample ceiling-exceedance share, which is zero for every feasible fit; see ceiling_exceedance(). print() of a CPB fit shows the median largest feasible count instead of the continuous ceiling, and describes alpha as the variance-to-mean ratio at an integer ceiling.
irr() takes its intervals as confint() does, reports odds ratios only under a logit link, leaves out fixed-effects dummies, and gives NA with a message for a term the data do not bound. For families whose exp(beta) is not a ratio of means it names the ratio ("lambda ratio", "centrality ratio", "long-run rate ratio", "location ratio", "mu ratio").
Stored bootstrap draws have one layout: named columns, the dispersion parameter on its natural scale (zi_gec() stored the log of delta; the draws of cpb_fe(), gec_fe() and zi_cpb() lacked it), and a fixed-effects fit’s reference-unit effect; saved objects keep their layout. The zero_fe dummies of zi_cpb() and zi_gec() are not bootstrapped and have no bootstrap standard error.
With cluster, the hurdles’ participation equation has cluster-robust standard errors (in hurdle_count(), those its se asks for).
Printed tables and irr() hide only the dummies of the fe, part_fe and zero_fe columns, so a factor() covariate in a formula is now shown. confint() of a cpb() fit names its columns "2.5 %" and "97.5 %", as stats::confint() does.
New refusals: first_difference() refuses a variable that also enters another model column (I(x^2), log(x), poly(), a spline), which gave a wrong contrast, and hurdle_gec() and zi_gec() refuse interactions; zi_test() refuses covariate-dependent inflation, a fixed-effects zi_cpb(), a separated zero equation, and a nest fitted to other rows, zero-truncated or on the other support; alpha_confint() refuses cpb_fe() fits, for which its interval was far too narrow.

New features

count_reg(), hurdle_count() and zi_count() gain the hyper-Poisson ("hyperpois", in the mean parameterization of Saez-Castillo and Conde-Sanchez 2013) and the exponentially weighted Poisson of Ridout and Besbeas (2004) with two or one tail weights ("ewp3", "ewp2"), with dhyperpois(), dewp() and their p/q/r functions, and the aliases "hyper_poisson", "ewp" and "exponentially_weighted_poisson". dispersion_test() tests the two weights of "ewp3" jointly (two-sided, by parametric bootstrap unless B = 0).
ame(): average marginal effects (a derivative or a contrast averaged over the rows of newdata) for every class, through the formula, so interactions and transformations are handled. Intervals follow first_difference(); B > 0 gives a refit-bootstrap interval, for a fixed-effects fit one for the average partial effect in the population of units. The derivative of a floor-support fit is refused, since its mean jumps at integer ceilings.
predict() gains se.fit, interval = "confidence" and level for the rows of newdata in every class (they were ignored), and offset = in every method. first_difference() gains at = to fix other model-matrix columns.
ceiling_exceedance() counts held-out observations above a CPB fit’s bound. predict.zi_cpb() gains type = "rate", and confint(fit, "alpha") gives a cpb() fit’s profile interval without a bootstrap.
The zero-inflated fits check their zero-inflation equation for separation, and the hurdles their participation equation, and warn when the data put no useful bound on a coefficient. The result is stored in separation (part_separation in the hurdles) and printed; a separated term has no standard error or interval. See ?zi_cpb, Separation.
compare_models() gains the columns truncated and separated, warns when it ranks a separated fit, a zero-truncated likelihood beside untruncated ones, or models fit to different responses, and notes when a single-equation fixed-effects fit sits beside a two-part fit that counts fewer unit effects.
tidy() gains scale (what each coefficient multiplies) and takes conf.int for every class. glance() gains BIC, dispersion, dispersion_parameter, dispersion_direction and coef_scale (and n.units for gec_fe()). Texreg columns are named by coefficient scale. augment() takes data =.
residuals(type = "pearson") for cpb_fe(), gec(), gec_fe() and count_reg(). count_reg() with robust or clustered errors also keeps the model-based ones, and summary() prints their ratio.
vcov() and confint() of a hurdle_cpb() or hurdle_gec() fit without intensity inference return the participation equation’s, with NA intensity rows, instead of NULL or an error.
A serial bootstrap that will take over a minute says so, fewer than 100 replicates bring a note on Monte Carlo error, and print() of dispersion_test() notes a p-value at its floor 1/(B + 1). compare_dispersion() gains converged, ud_screen() names comparators it did not fit, and summary() of cpb() and count_reg() fits counts the rows dropped for missing values.

Improvements

The zero-inflated fits start from more points (such as the zero-inflated Poisson’s coefficients and a low-inflation start), followed by a BFGS pass, so they can reach higher maxima. zi_cpb(fe =, method = "em") adds block coordinate ascent on the observed-data likelihood, and gec() and King’s cpb() also end with a BFGS pass.
count_reg(), the intensity of hurdle_count() and the starts of zi_count() evaluate the likelihood over distinct rows weighted by their counts: the same results, 5 to 20 times faster with repeated rows.
Every class bootstraps through one engine (the same resampling, failed or non-converged replicates dropped with a warning, a warning for few clusters). Parallel workers load the package from the session’s library.
More argument checks (missing values in time, vectors of the wrong length, offset in score() without newdata). count_reg(family = poisson) accepts the family object, and cpb_fe() and gec_fe() say where two-way fixed effects go.
Documentation: each family’s coefficient scale (?irr.count), sources and limits; the robust Poisson as the benchmark for mean effects, and when the CPB’s are consistent (?count_reg, ?cpb, a new vignette section); what the implied ceiling bounds and why the held-out check is one-sided (?implied_ceiling, ?ceiling_exceedance); King’s notation and why Hessian standard errors fail (?cpb); incidental-parameters bias in fixed-effects slopes and participation logits (?cpb_fe, ?hurdle_cpb); what dispersion_profile() can and cannot show.

Bug fixes

cpb_fe() and gec_fe() did not reorder a vector cluster with the rows, which put the wrong clusters on panels not supplied unit by unit. A cluster factor with unused levels no longer draws empty clusters.
A fixed-effects unit with no positive count sits 30 log units below its own largest linear predictor, not at an absolute -30 (gec_fe(), cpb_fe() on King’s support, zi_cpb(fe =), dispersion_test()’s null); with an uncentered covariate the absolute value gave NA fitted values and shifted the estimates and the test statistic.
gec_fe() refuses a likelihood term that the support guard cuts, as gec() does, instead of renormalizing, which overstated the likelihood.
The reference unit of a fixed-effects fit, used in first differences and for unseen units, averages only the units whose effect the data identify (leaving out units with no positive count and, in a zero-truncated King’s CPB, units whose counts are all 1); these units had moved it by up to an order of magnitude. Dummy-coded unit effects follow the same rule.
An unseen unit, or newdata without the unit column, enters the part_fe and zero_fe equations and the count equations of count_reg(fe =), hurdle_count(fe =) and zi_count(fe =) at the reference unit in predict(), score() and cv_score(), and held-out scoring treats an unseen level of any factor the same way. Several classes used to stop there.
In new data the two-part fits evaluate a unit at a limit at that limit: in the hurdles a unit that never participates gets participation probability zero and one that always participates gets one; in the zero-inflated fits a zero_fe unit gets structural-zero probability one if it has no positive count and zero if every count is positive, and a count-equation unit with no positive count a mean of zero. They used the optimizer’s finite stopping values, which a covariate could move to almost any probability.
zi_cpb(fe =) stops when its first-stage fit fails (it continued with every unit effect at zero), and its warning that the count component cannot cover the data no longer fires on valid fits.
Delta-method intervals step each parameter by its standard error, not by a share of its value, so they do not depend on covariate units (a covariate in small units gave intervals of 1e80 or NA).
first_difference(B > 0) on hurdle_cpb() and zi_cpb() refits each replicate with the fit’s own settings and resamples its units or clusters (it refit the formulas alone); an offset or weights must then be columns of data. On hurdle_cpb() it uses the fitted link, not always the logit, and the two-part classes hold an offset at its mean.
count_reg(fe =), hurdle_count(fe =) and zi_count() keep a transformed response and accept y ~ .; a | inside I() is no longer read as a two-part formula; character covariates predict on new data that hold only some of their values.
The gamma count’s log-probabilities no longer underflow at large means, its moments are refused at the support cap, and its variance no longer overflows. The variances of the COM-Poisson, double Poisson and generalized Poisson are computed in two passes, which fixes large errors for concentrated distributions.
The COM-Poisson (both parameterizations), double Poisson and generalized Poisson below zero return NA where the 1e5-term cap would truncate their normalizing sum (a mean of order 1e5), and the likelihood treats such points as infeasible; they returned an overstated likelihood. This also stops very slow mean-parameterized COM-Poisson fits.
The mean-parameterized COM-Poisson finds its rate at a mean below 1 and a large nu, and its root search for nu < 1 starts between the Poisson and geometric rates. The COM-Poisson’s nu and the hyper-Poisson’s gamma are held in [1e-4, 1e4]; a fit at that edge is compared with one at the family’s limit, and print() notes it.
The generalized Poisson with lambda >= 0 uses its closed-form normalizer (the capped sum overstated the log-likelihood for long tails), and its predicted probabilities no longer fail for a slightly negative lambda at large means.
rgec() and qgec() return NA with a warning where the support guard binds, a non-finite rate gives NA GEC moments, and the truncated means of hurdle_gec() and hurdle_count() share one floor on 1 - P(0).
Failures are reported, not masked: a covariance that cannot be formed, or a singular information, warns; a non-finite gradient of the concentrated fixed-effects objective falls back to central differences (zeros stopped BFGS); the matched families report convergence only when the BFGS stage or its polish converged or the fit is saturated; a texreg registration failure is reported at load.
glance(), the BIC of compare_models() and texreg’s “Num. obs.” use a weighted fit’s weight total; max.support is capped one below the integer maximum; the help pages no longer lose text to escaped percent signs.

underdisp 0.1.2

CRAN release: 2026-09-30

Bug fixes

Converting a fitted ceiling to an integer is now guarded against values outside the range of int. The CPB’s ceiling is lambda / (1 - alpha) and the GEC’s is mu / (1 - delta); both diverge as the dispersion parameter approaches its equidispersed limit, which an optimizer reaches routinely, and the conversion was undefined behaviour there. CRAN’s macOS sanitizer build reported it in cpb_pmf.h and cpb_fe.cpp. All six conversions of this kind, including two the sanitizer did not reach, now cap the ceiling at max.support + 1 first, through a single rule in src/support_cap.h.

A ceiling past max.support was already treated as infeasible, on both of the platforms whose out-of-range conversions differ, so fitted values are unchanged: every cross-sectional and fixed-effects CPB fit, every GEC fit without fixed effects, and the calibrated interval reproduce bit for bit. The exception is gec_fe() on data whose dispersion is driven to the Katz limit, where the breakpoint scan now runs from a defined bound and the optimizer can settle a few parts in 10^8 away.

underdisp 0.1.1

CRAN release: 2026-09-28

New families

count_reg(), hurdle_count(), and zi_count() gain four underdispersed families through the same interface, each inheriting fixed effects, zero-truncation, hurdle and zero-inflated forms, offsets, weights, the sandwiches, proper scores, DHARMa simulation, and the QoI methods: "mpcmp" (the COM-Poisson in Huang’s 2017 mean parameterization, with the rate solving E(Y) = mu found in C++), "genpois" (the Consul-Jain generalized Poisson with constant dispersion lambda, renormalized on its finite support when lambda < 0), "gammacount" (Winkelmann’s 1995 renewal-process model with gamma waiting times), and "doublepois" (Efron’s 1986 double Poisson with the exact normalizing constant). Every family’s pmf was checked against an independent implementation (rmutil, gamlss.dist, VGAM) to machine precision, and the generalized-Poisson and mean-parameterized COM-Poisson log-likelihoods reproduce glmmTMB’s genpois and compois fits exactly.
dgammacount(), ddoublepois(), dgenpois() and their p/q/r companions; the COM-Poisson functions accept mu = for the mean parameterization.

New diagnostics

dispersion_test(): the likelihood-ratio test of a fitted family’s dispersion parameter against its Poisson value on the fitted design (fixed effects included, the unit effects profiled out of the Poisson null), with two-sided or directional alternatives by the signed root and the boundary mixture where the Poisson value is on the boundary (CPB, negative binomial); and the Cameron-Trivedi auxiliary-regression test for a Poisson fit. The p-value is a parametric bootstrap under the fitted Poisson (B, cores) whenever the Poisson value is on the boundary of the family’s parameter space (CPB, negative binomial) or the fit has unit fixed effects; where it is interior, the asymptotic p-value is used unless the first-order bias from the null’s estimated mean parameters, p / sqrt(2n) standard deviations, pushes a 5% test’s size above 6%. summary() of a CPB fit flags its asymptotic p-value.
ud_screen(): in ztp_threshold = "bootstrap" mode the NB-vs-Poisson p-value is a parametric bootstrap under the fitted Poisson (nb_boot) when the mean model is within the comparator gates; the fast mode’s asymptotic value is labeled conservative. The bootstrap restores the random-number state, so the at-risk threshold and every other result are unchanged.
dispersion_profile(): the empirical conditional variance-to-mean ratio by bins of the fitted mean against the ratio each fitted family implies, with a plot() method, so the mechanism behind an underdispersed outcome (a hard ceiling, regular event timing, a soft tail) can be read off the data.
calibrate_alpha(): an opt-in calibration of the interval for the CPB dispersion parameter. The support of the CPB moves with the parameters through the log-ceiling coefficients, which behave like endpoint parameters, and the estimate of alpha inherits a bias toward zero from them; the first-order (chi-square) profile interval therefore covers about 0.90 to 0.94 in simulations from the CPB, with its misses on the upper side. calibrate_alpha(fit, B, cores, seed) stores the parametric-bootstrap distribution of the signed root of the profile likelihood ratio in the fit, and alpha_confint(), confint(), implied_ceiling() and summary() then report the interval cut at its quantiles (0.92 to 0.96 in the same simulations, misses balanced). The responses are drawn once on the calling process and the refits use no random numbers, so the result does not depend on cores. It is refused, with the reason, for a fit at the max.support guard, weights that are not whole numbers, and a fit saved by 0.1.0. With few distinct covariate patterns the first-order interval already covers at its level and the calibration is not needed. Without it the accessors report the first-order interval, and attr(, "method") and summary() say which one is shown.

Weights, parallel bootstraps, and column scaling

weights = (frequency weights) on every estimator; a weight of w is equivalent to w copies of the row in the likelihood, the information, the score products, and the bootstrap, and nobs() returns the weight total. In the parametric bootstraps (dispersion_test(), calibrate_alpha()) a row of weight w contributes w independent simulated responses; weights that are not whole numbers have no such replicate data, so dispersion_test() then reports the asymptotic p-value with a note.
cores = on every bootstrap, the parametric-bootstrap screen threshold, cv_score(), and the refit first differences: replicates run on a socket cluster seeded from the calling session.
Every optimizer works on covariates scaled to unit standard deviation and maps the coefficients back, and the Nelder-Mead fits restart from their own solution until the improvement stops, so a fit no longer depends on the units of a covariate. The numerical Hessians and score matrices of the count families are differenced in the same scaled parameterization.

The CPB mean

fitted(), predict(type = "response"), residuals(), and first_difference(quantity = "mean") for cpb, cpb_fe, and the hurdle-CPB intensity now report the exact mean of the fitted distribution (a finite sum of the renormalized pmf, and the conditional mean E(Y | Y >= 1) for zero-truncated fits), as gec() and the matched families already did. The rate parameter exp(x’b) equals the mean only when lambda/(1-alpha) is an integer; it remains available as $rate and predict(type = "rate"). The same holds for the generalized event count: delta is the Katz dispersion parameter, which is the variance-to-mean ratio on an unbounded support, and exp(x'b) is the rate of the recursion.
Zero-truncated count_reg() fits report E(Y | Y >= 1).

Methods

confint() for every model class (bootstrap percentile, normal-approximation from the bootstrap, or Wald from the analytic covariance, one rule per inference type); vcov() for the hurdle and zero-inflated classes; implied_ceiling() for cpb_fe; irr() and first_difference() on a cpb fit without a bootstrap return point estimates with method = "none"; vcov() returns NULL rather than an error when no inference was requested.
irr() labels the binary stage "participation" (hurdle) or "inflation" (zero-inflated), which point in opposite directions.
texreg extract() methods for count_reg, hurdle_count, zi_count; the methods for every class now attach to texreg’s generic (in 0.1.0 the registration silently failed and tables fell back to broom), and the count families label their dispersion row by the parameter.
tidy() for the two-part classes prefixes each term by its component, so modelsummary() lays the table out without a shape argument.
predict(), score(), and the QoI methods rebuild a newdata model matrix in the fit’s own factor coding (stored levels and contrasts) and align it to the coefficients by name; a factor whose reference level was set through its contrasts attribute (as in AER’s NMES1988) is no longer re-coded on newdata, and a coding mismatch is refused instead of multiplied by position.
A formula offset() term and a pscl-style y ~ x | z formula are refused with a message naming offset =, participation =, and zero =; rows whose unit identifier is missing are dropped like other incomplete rows; the screen’s generalized-Poisson comparator is the package’s own family, so it can report underdispersion, and it is gated by comp_max_par/comp_max_n like the COM-Poisson (a full fit, not a diagnostic, on a dummy-heavy design); print.gec() labels the dispersion direction by the likelihood-ratio test; count_reg(fe =) hides the absorbed dummies and warns when cluster is the fixed-effects variable; standard errors from a Hessian that is not positive definite are NA with a warning; the generalized Poisson’s support is capped for negative dispersion; score() honours frequency weights; the auxiliary test is weighted; the parallel driver draws its cluster seed eagerly, so consecutive cores > 1 bootstraps differ.
converged is stored by every class and reported by print()/summary() when it is FALSE; a matched-family fit whose optimizer never left the infeasible plateau (a fitted mean outside [1e-12, 1e8]) is reported as not converged with a warning, and its log-likelihood is -Inf, instead of returning the sentinel as a likelihood.
The predictive distribution of a zero-truncated gec() or count_reg() fit, as score(), rootogram(), pit_hist(), and cv_score() use it, is the fit’s own zero-truncated pmf normalized over the full support (in 0.1.0 the Katz pmf was scored untruncated and the count families renormalized over the observed range), so the in-sample log score equals -logLik/n for every class; fitted() on a hurdle_gec is the exact conditional mean that predict() returns.
cv_score(), rootogram(), pit_hist(), and dispersion_profile() apply frequency weights as score() does, so a row of weight w counts as w copies of itself: the rootogram’s observed and expected frequencies, the PIT histogram, the dispersion profile’s bins (equal weight, rows kept whole) and bin means, and, through a new weights argument, the mean held-out scores of cv_score() and score(newdata =). rootogram() no longer fails when kmax exceeds the largest count.
pit_hist() read the lower bound F(y - 1) out of row order wherever a zero count preceded a positive one, so its histogram was wrong for data with zeros; the bound is now looked up row by row.
ud_screen() reports the NB-vs-Poisson likelihood ratio at its boundary value 0 when glm.nb() stops short of the Poisson limit, and prints “at the Poisson boundary” there rather than an overdispersion reading.
The zero-truncated CPB at an implied ceiling below 1 is the point mass at 1 in the pmf, the moments, the first differences, and the simulator alike (the zero-truncated distribution for every ceiling in [1, 2)), so the hurdle’s predictions, scores, and rootograms are finite at rows whose covariates put the rate below a unit’s feasibility floor.
One rule names the coefficients of the two-part classes everywhere: participation:/intensity: for the hurdles and count:/zero: for the zero-inflated models, in vcov(), confint(), tidy(), and the texreg tables, so a term read off a table indexes the covariance and the interval.
A vector cluster is checked for its length and for missing values by every estimator before it is used (the fixed-effects and two-part estimators recycled or padded a short vector in 0.1.0), and a row dropped for a missing unit identifier drops its entry of a vector offset, weights, or cluster with it; hurdle_count() selects the cluster-robust covariance when cluster is supplied, as its siblings do, and accepts a vector cluster aligned to the data; zi_cpb() and zi_gec() refuse a non-integer response instead of flooring it; zi_cpb(), zi_count(), and zi_gec() work on scaled covariates with a Nelder-Mead polish like every other estimator (the unscaled joint fit of zi_count() could stop at its starting values: on lme4’s grouseticks it stopped 4.4 log-likelihood units below pscl’s optimum); summary.zi_count() prints the stored inflation standard errors; every bootstrapping summary states how many replicates converged; count_reg() and zi_count() store $fitted.values on the mean scale and the natural parameter as $mu.
The parallel workers of cores > 1 attach the package, so a cv_score() fit function that names cpb() or count_reg() unqualified runs on them.
The analytic, robust, and cluster covariances of the matched families are differenced at an absolute step of 1e-3 in the scaled parameterization (numDeriv’s default relative step left the family’s support at large fitted means and returned NA standard errors).
predict() with newdata returns NA for a row whose covariate or offset is missing, as predict.glm() does, instead of failing inside the distribution routines. A column holding only missing values (which R reads as logical) is treated as a missing value of the fitted variable’s type, and a variable supplied with a different type than in the fitting data is refused with a message that says so.

Estimation and numerics

The CPB profile-likelihood interval for alpha is traced by continuation with a feasibility-repaired start at each alpha, and the profile at each alpha is maximized with the fit’s own multistart (BFGS, Nelder-Mead polish, perturbed restarts), so the lower limit is the re-maximized profile’s crossing point; the interval is one-sided only when the profile has not dropped to the cut by alpha = 0.005. Each limit is the outermost crossing of its cut (the profile is a saw-tooth, so the set above a cut need not be connected), and a cut deeper than the stored trace, as at a level above 0.995, re-runs the trace instead of returning the parameter bounds.
rcpb(truncated = TRUE) warns when a rate’s ceiling lambda / (1 - alpha) is below 1: the zero-truncated CPB does not exist there, and the draw of 1 it returns has probability zero under the stated parameters.
max.support defaults to max(500, 10 * max(y)); a fit whose largest implied ceiling reaches the guard warns and records $support_binding; a fit whose every start is infeasible errors instead of returning the starting values; a CPB fit that the guard stops short of its Poisson limit reports the likelihood-ratio statistic at its boundary value 0.
The concentrated fixed-effects fits solve each unit’s intercept exactly. The unit objective is a saw-tooth in the intercept (the support floor moves as the rate crosses integer multiples of 1 - alpha, and the objective drops at every such breakpoint), so its supremum sits at an interior critical point of a tooth or at the left limit of a breakpoint; the inner search enumerates the breakpoints within a window of the envelope’s peak (as wide as the widest tooth in the unit) and refines inside whole teeth by golden section; every evaluation searches each unit afresh, so the concentrated likelihood is a function of the parameters alone. It is continuous in the dispersion parameter but can jump in the slopes where a unit sits on its feasibility floor, so the outer optimizer is a deterministic multistart: BFGS with the analytic envelope gradient (the breakpoint a unit sits at is tracked in the gradient), a Nelder-Mead polish, and two perturbed restarts; ?cpb_fe states the size of the alternative optima. The Katz family below delta = 1 has a finite support too, but its likelihood only kinks at an integer ceiling (the entering point’s mass grows from zero), so gec_fe() maximizes each unit’s unimodal objective by golden section on an expanding bracket, with the feasibility floor exact (support 0..ceiling(mu/(1-delta))). The compiled code asserts the unit index against the design, and the fits use only the estimation rows, so rows dropped for missing values no longer misalign the unit blocks.
cpb() documents the discontinuity of the CPB log-likelihood and the deterministic multistart that maximizes it (BFGS from the Poisson solution at each alpha.start, Nelder-Mead polish, feasibility-repaired restarts); the mean-parameterized COM-Poisson tabulates its log-factorials.
The CPB pmf is normalized by a recursion outward from the mode of its terms (consecutive terms differ by the factor (n - k)/(k + 1) (1 - alpha)/alpha, so the sum runs over the body of the distribution and stops where the terms fall below relative precision on either side), shared by the pooled likelihoods and the concentrated fixed-effects likelihood and gradient. A likelihood evaluation costs the width of the pmf rather than its ceiling: on counts in the hundreds one evaluation of the concentrated objective takes 0.04 s instead of 2.6 s, and the pooled objective is 45 times cheaper. The R implementation of the pmf keeps the direct lgamma form as the independent check (agreement to 5e-10 across alpha from 0.01 to 0.999 and rates from 0.05 to 2,000).
The reported maximum of a pooled cpb() fit is the maximum of the fit’s own profile in alpha: the profile is traced by continuation on both sides of the estimate (the slopes re-maximized at each step from the neighbouring solution and from the Poisson slopes, with a finer pass next to the estimate) until it has dropped four log-likelihood units, and a trace point above the multistart’s value restarts the fit from there; the multistart itself runs from nine values of alpha and ends with an exact breakpoint search in alpha at the fitted slopes, since the teeth of the surface in alpha are thinner than any grid. confint() and summary() read the profile interval off the stored trace. On small samples the multistart alone could stop on a lower tooth with alpha off by 0.1; the remaining differences against a dense grid of starts are of the size of the jumps (below 0.005 log-likelihood units on 12 samples of 72 zero-truncated counts). Unit intercepts belong in cpb_fe(), which solves them exactly; a pooled fit with many dummy columns can stop short of it (0.14 log-likelihood units on pscl’s prussian with 14 corps dummies).
The inner search of cpb_fe() locates the peak of a unit’s smooth part to 1e-6, so the right tooth is searched at ceilings above a few hundred, and searches the neighboring tooth when the interior optimum sits at a tooth’s end.
The GEC support guard refuses instead of renormalizing: a rate whose recursion has not reached the tail by max.support is infeasible in the likelihood, a fit whose fitted support reaches the guard warns and records $support_binding, and predict(), dgec(), and pgec() return NA with a warning for such a rate (in 0.1.0 the pmf was silently renormalized over 0..max.support, so a prediction far outside the fitted range was too small).
The COM-Poisson support rule no longer clips the rate at 1e10 (with mu = at strong underdispersion the rate exceeds it and the pmf was wrong), and the d/p/q/r functions return NA for a non-finite rate.
The final Nelder-Mead polish chains its restarts until a restart gains less than the tolerance, as documented.
Design matrices are rank-checked (an empty factor level or an exact linear combination is refused with the column named; a covariate constant within units is refused by the fixed-effects fits); a wrong-length vector offset, weights, or cluster is refused; cluster supplied without a bootstrap is announced; a cluster bootstrap with one cluster is refused and with fewer than thirty warned.
score() raises kmax to the largest held-out count, and its predicted distributions carry the fit’s offset (the in-sample log score equals -logLik/n for every class); the hurdle intensities are zero-truncated over the full support; first_difference() holds an offset at its mean and refuses a variable that enters an interaction; the Katz recursion runs at least to the observed count before its tail is dropped, so far-tail counts under an unbounded member get their small probability rather than an infeasibility sentinel; integer ceilings are feasible (floor with a tolerance) in the C++ and the R pmf alike; the cluster-robust sandwich carries the G/(G-1) factor.
The two-part models align both equations on the complete cases of every model variable; zi_cpb() records the user’s call; predict() on gec and gec_fe uses the stored factor levels; the exported d/p/q/r functions return numeric(0) for zero-length input and refuse a vector dispersion parameter.
Every estimator drops unused factor levels from its model frame, as lm() and glm() do, so data subset from a larger set no longer trip the rank-deficiency refusal on an empty level.
A response whose counts are all zero still fits, with a warning that the rate runs to its lower bound and the dispersion parameter is not identified.
max.support, inner_it, alpha.start, and B (with se = "bootstrap") are checked on entry and refused by name before any fitting: a support below 1 or beyond the integer range, a fractional or negative count, fewer than two bootstrap replicates, a start outside (0, 1).
The support guard’s warning comes once, from the fit the user asked for: the refits of a bootstrap, a bootstrap first difference, or the split-panel jackknife, the start values of the zero-inflated fits, the nest that zi_test() fits from a formula, and the hurdle arm of compare_dispersion() no longer repeat it, and a gec() or gec_fe() fit no longer issues it twice.

Package

weights, cores, and unknown-argument checks on the matched families (misspelled arguments are errors); inst/CITATION; a pkgdown site at https://bagozzib.github.io/underdisp/ and GitHub Actions running the full test suite (including the batteries skipped on CRAN) on five platforms.

underdisp 0.1.0

CRAN release: 2026-08-20

First public release. underdisp provides diagnostics and a unified family of estimators for underdispersed count data (conditional variance below the conditional mean), the case the Poisson and negative binomial defaults cannot represent. The likelihood core is implemented in C++.

Screening

ud_screen(): at-risk dispersion diagnostic that benchmarks a fitted count against a zero-truncated Poisson under unit and period fixed effects, with a marginal statistic, a zero-truncated-Poisson at-risk statistic, and an over-conditioning guard. The at-risk threshold is available in two modes, ztp_threshold = "calibrated" (a fast n-dependent calibrated cutoff) and ztp_threshold = "bootstrap" (a parametric bootstrap of the fitted zero-truncated-Poisson null on the data’s own design, ztp_boot_B), and an optional native COM-Poisson comparator (run_comp).

Estimators

cpb(): continuous parameter binomial (CPB) regression and its zero-truncated variant, with an interpretable observation-specific bound (the implied ceiling). C++ multistart estimation.
cpb_fe(): high-dimensional fixed-effects CPB via a concentrated likelihood that matches the full-dummy likelihood to machine precision.
hurdle_cpb() and zi_cpb(): two-part (participation logit + zero-truncated CPB intensity) and true zero-inflated mixture forms, each with fixed effects available in both equations.
gec(): King’s generalized event count / Katz-family regression, a single model whose dispersion is estimated freely and spans underdispersion, equidispersion, and overdispersion, with gec_fe(), hurdle_gec(), and zi_gec() variants.
count_reg(), hurdle_count(), zi_count(): matched Poisson, negative binomial, and COM-Poisson families sharing one interface, with fixed effects (including two-way), offsets, and truncated/hurdle/zero-inflated forms, for like-for-like comparison against the CPB and GEC families.

Inference

Coefficient inference by a cold-multistart pairs (or cluster) bootstrap with percentile intervals; dispersion inference by a profile-likelihood interval for the shape parameter, which respects the feasibility boundary where the bootstrap distribution is skewed (alpha_confint()).
Robust and clustered standard errors throughout, via a block bootstrap rather than a Hessian sandwich, since the CPB’s parameter-dependent support makes the Hessian unreliable. Fixed-effects bootstraps use normal-approximation intervals, because replicates share the incidental-parameters bias and percentile intervals would be off-center.

Bias-corrected fixed-effects estimation

bias_correct = "jackknife" on the concentrated fixed-effects estimators (cpb_fe(), gec_fe(), and the fe = intensity of hurdle_cpb()): the split-panel jackknife (Dhaene & Jochmans 2015) removes the leading 1/T incidental-parameters bias of the dispersion parameter, validated by the package’s fixed-effects bias Monte Carlo (alpha bias at T = 6/10/20/40 falls from -.147/-.096/-.056/-.033 to -.038/-.022/-.016/-.009). The correction carries a validity gate for the method’s time-homogeneity requirement: the time-split disagreement is compared against exchangeable unit-split placebo noise, and on panels where the two halves do not estimate a common parameter the correction is refused with an informative warning and the maximum-likelihood fit is returned. A refusal is diagnostic information about the panel, not an error.
For separation in dummy-heavy participation equations, the hurdle factorizes: fit the participation stage with a dedicated bias-reduction package (logistf, brglm2) alongside the zero-truncated intensity from this package.

Comparison and testing

compare_dispersion() and compare_models(): information-criteria and proper scoring-rule comparison across Poisson, negative binomial, COM-Poisson, and CPB (and their hurdle and zero-inflated forms), with type-safe rejection of non-package or mismatched objects.
zi_test(): boundary-corrected likelihood-ratio test for zero-inflation (the plain model nested in the mixture at the boundary), for both the CPB and GEC families.

Quantities of interest

predict(), implied_ceiling(), alpha_confint(), irr(), and first_difference(). Every first_difference() method in the package returns one unified "ud_fd" data frame (component, from, to, diff, lower, upper, method): single-equation fits report the mean, and the two-part/mixture models report the exact extensive/intensive decomposition (binary-stage probability, count-stage mean, marginal mean), with a stage argument on every two-part method to move the covariate in a single equation and a per-row method label recording the uncertainty source (stored bootstrap for cpb, refit bootstrap via B= for the CPB two-parts, delta-method where an analytic covariance exists). Unknown arguments error rather than being silently ignored.

Distributions and simulators

Density, distribution, quantile, and random-generation functions for the package’s own families: dcpb/pcpb/qcpb/rcpb, dgec/pgec/qgec/rgec, dcompois/pcompois/qcompois/rcompois, plus the mixture simulators rhurdle_cpb() and rzicpb(). The standard zero-truncated, hurdle, and zero-inflated Poisson/negative-binomial forms are deliberately internal: CRAN already serves those families (e.g. VGAM), and the exported surface is reserved for what the package adds.
simulate() methods for every model class (cpb, cpb_fe, gec, gec_fe, count_reg, hurdle_count, zi_count, hurdle_cpb, hurdle_gec, zi_cpb, zi_gec), with fitted() completed across the two-part and fixed-effects classes, so any fit plugs into DHARMa’s simulated-residual diagnostics via DHARMa::createDHARMa().

Panel tools and calibration

mundlak(): correlated-random-effects device that augments a formula with unit means, model-agnostic across the estimators.
score() and cv_score() (log score and ranked probability score, in sample and by cross-validation), rootogram() (Tukey hanging), and pit_hist() (non-randomized PIT histogram).

Numerical robustness

The estimators are hardened for the difficult likelihoods this package targets (near-boundary dispersion, zero-heavy and dummy-heavy panels):

Multistart estimation with per-start restart-polishing, and, for the zero-inflated CPB, an expectation-maximization cross-check run alongside the direct maximum-likelihood multistart so the better optimum is kept, avoiding the degenerate zero-mass-collapse basin.
Loud, estimator-backed fallbacks that degrade gracefully rather than returning NA: a Fisher-scoring fallback for the zero-truncated-Poisson screen when the default IRLS fit errors or silently diverges, and a ridge-regularized Fisher step in the bootstrap null on singular information matrices. Each fallback warns and reports the method actually used.
Fitted-rate sanity guards that key on the pathology (diverged rates) rather than its correlates (iteration counts), and floored mixture likelihoods so a transient zero-probability iterate cannot produce a non-finite objective.
Fit-rank-exact degrees of freedom in the at-risk screen and the mixture comparisons, so zero-heavy designs are neither over- nor under-penalized.

Data and integration

Bundled peacekeeping dataset (UN peacekeeping contributions, a state-year roster with structural zeros).
broom tidy()/glance() methods for all model classes (augment() for cpb), and modelsummary / texreg table support.