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.
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.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 = (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.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.count_reg() fits report E(Y | Y >= 1).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.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.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.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.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.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.cores > 1 attach the package, so a cv_score()
fit function that names cpb() or count_reg() unqualified runs on them.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.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.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.(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).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).cpb_fe() locates the peak of a unit's smooth part to
1e-6 (in 0.1.1's first candidate the twelve golden-section steps left an
error of 1.6e-3, wider than a tooth at ceilings above a few hundred, so the
wrong tooth could be searched) and searches the neighbouring tooth when the
interior optimum sits at a tooth's end.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).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..ud_nm_polish() chains its restarts until a restart gains less than the
tolerance, as documented (the first candidate stopped after one).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.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.lm()
and glm() do, so data subset from a larger set no longer trip the
rank-deficiency refusal on an empty level.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).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.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.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++.
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).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.alpha_confint()).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.logistf, brglm2) alongside the zero-truncated intensity from
this package.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.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.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().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).The estimators are hardened for the difficult likelihoods this package targets (near-boundary dispersion, zero-heavy and dummy-heavy panels):
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.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.