underdisp provides tools for detecting and modeling
underdispersion in count data: the case where the conditional
variance is below the conditional mean, so counts cluster more tightly
around their expectation than a Poisson allows. The Poisson and negative
binomial defaults cannot represent it; the negative binomial in
particular collapses onto the Poisson when the data are
underdispersed.
We generate a count with a conditional variance-to-mean ratio of about one half.
ud_screen() returns a marginal verdict, and, for
zero-inflated outcomes, an at-risk verdict benchmarked against a
zero-truncated Poisson (which is what separates genuine
underdispersion from the artifact of conditioning on positive
counts).
ud_screen(y ~ x, data = d, run_cpb = FALSE)
#>
#> === Underdispersion screen ===
#> Formula: y ~ x
#> N=400 mean=5.527 max=22 %zero=1.8 (unconditional var/mean=1.89)
#>
#> MARGINAL verdict: UNDERDISPERSED
#> Pearson=0.564 prop.slope=-0.391 (p=<2e-16)
#> NB vs Poisson LR = -0.01 (p= 0.5 ; sig => overdispersion)
#>
#> AT-RISK (y>0) verdict:UNDERDISPERSED [n_pos=393]
#> ZTP-Pearson = 0.569 (underdispersed if < 0.885, the calibrated 5% threshold)
#>
#> Model comparison (log-lik):
#> Poisson NB GP COMP CPB
#> -798.901 -798.904 -798.901 -772.010 NA
#>
#> COM-Poisson (full data): nu = 1.84 (> 1 = underdispersed, soft tail)cpb() fits the CPB, with truncated = TRUE
for the common case in which underdispersion lives among the positive
counts of a zero-inflated outcome.
fit <- cpb(y ~ x, data = d[d$y > 0, ], se = "none")
summary(fit)
#>
#> Continuous Parameter Binomial regression (zero-truncated)
#> N = 393 inference: none
#>
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.5797 NA NA NA
#> x 0.4931 NA NA NA
#>
#> alpha = 0.5279 (profile 95% CI: 0.509 to 0.613)
#> Implied ceiling lambda/(1-alpha): median 10.24 range 2.47 to 37.96
#> logLik = -745.06 AIC = 1496.13
#> LR vs ZT-Poisson (H0: alpha = 1): 58.62, p 9.5445e-15The dispersion parameter alpha summarizes the
compression, and each observation carries an implied ceiling
lambda / (1 - alpha).
Predicted probabilities, the implied ceiling, and King-style first differences are available for user-specified covariate profiles.
Because the CPB’s support depends on its parameters, Hessian-based standard errors are unreliable; coefficient inference uses a cold-multistart pairs bootstrap (validated to nominal coverage in the companion paper), and the dispersion parameter carries a profile-likelihood interval.
fit_b <- cpb(y ~ x, data = d[d$y > 0, ], se = "bootstrap", B = 99)
summary(fit_b)
#>
#> Continuous Parameter Binomial regression (zero-truncated)
#> N = 393 inference: bootstrap
#>
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.579699 0.017471 90.417 < 2.2e-16 ***
#> x 0.493096 0.018320 26.915 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> alpha = 0.5279 (profile 95% CI: 0.509 to 0.613)
#> Implied ceiling lambda/(1-alpha): median 10.24 range 2.47 to 37.96
#> logLik = -745.06 AIC = 1496.13
#> LR vs ZT-Poisson (H0: alpha = 1): 58.62, p 9.5445e-15
#> (99 bootstrap resamples converged)
irr(fit_b) # incidence-rate ratios with percentile intervals
#> term equation ratio estimate lower upper method
#> (Intercept) count IRR 4.853 4.722 5.045 bootstrap (stored)
#> x count IRR 1.637 1.579 1.688 bootstrap (stored)
alpha_confint(fit_b) # profile-likelihood interval for alpha
#> lower upper
#> 0.5093570 0.6132917
#> attr(,"alpha")
#>
#> 0.5278969
#> attr(,"boundary")
#> boundary
#> FALSEgec() fits King’s generalized event count (Katz) model,
whose dispersion delta (the variance-to-mean ratio) is
estimated freely — so the data choose the direction of dispersion rather
than the analyst presuming it. On the underdispersed count above it
recovers delta well below one; on a Poisson outcome it sits
at one.
gec(y ~ x, data = d, se = "none") # delta ~ 0.5
#> Generalized event count (Katz family) regression
#> Call: gec(formula = y ~ x, data = d, se = "none")
#> (Intercept) x
#> 1.5768 0.4964
#>
#> dispersion delta (Var/Mean) = 0.553 [underdispersed]
#> logLik = -769.78, n = 400
gec(y ~ x, data = data.frame(y = rpois(n, exp(1 + 0.4 * x)), x = x),
se = "none") # delta ~ 1
#> Generalized event count (Katz family) regression
#> Call: gec(formula = y ~ x, data = data.frame(y = rpois(n, exp(1 + 0.4 * x)), x = x), se = "none")
#> (Intercept) x
#> 0.9650 0.3493
#>
#> dispersion delta (Var/Mean) = 0.954 [underdispersed]
#> logLik = -740.83, n = 400The GEC carries the same zero-truncated, hurdle
(hurdle_gec()), zero-inflated (zi_gec()), and
fixed-effects (gec_fe()) variants as the CPB.
Underdispersion is typically a within-unit phenomenon that
pooled analyses hide. cpb_fe() absorbs a full set of unit
fixed effects by concentrating them out of the likelihood, so it scales
to thousands of units.
panel <- do.call(rbind, lapply(1:50, function(i) {
xx <- rnorm(12); NN <- pmax(round(exp(rnorm(1, 0, 0.4) + 0.4 * xx) / 0.5), 1)
data.frame(unit = i, x = xx, y = rbinom(12, NN, 0.5))
}))
cpb_fe(y ~ x, data = panel, fe = "unit")
#> CPB regression with 50 unit fixed effects (concentrated likelihood)
#> Coefficients:
#> x
#> 0.3119
#>
#> alpha (shape parameter): 0.4769 median implied bound: 2.27
#> Note: alpha is subject to incidental-parameters bias for short panels; see ?cpb_fe.compare_dispersion() fits the Poisson, negative
binomial, the native soft-tail COM-Poisson, the free-dispersion GEC, and
the hard-ceiling CPB, and reports a fit comparison plus the CPB’s
ceiling-exceedance share.
compare_dispersion(y ~ x, data = d)$table
#> df logLik AIC BIC logscore rps zero_fit
#> Poisson 2 -798.9009 1601.802 1609.785 1.997252 0.9855061 0.024243957
#> NegBinomial 3 -798.9042 1603.808 1615.783 1.997261 0.9855095 0.024244780
#> COM-Poisson 3 -772.0096 1550.019 1561.994 1.930024 0.9631112 0.009312096
#> CPB 3 -769.2026 1544.405 1556.380 1.923006 0.9642624 0.010309391
#> GEC 3 -769.7769 1545.554 1557.528 1.924442 0.9644021 0.010538704For model selection, count_reg() fits Poisson,
negative-binomial, and COM-Poisson regressions – each with the same
fixed-effects, zero-truncation, hurdle, zero-inflation, offset, and
robust-/cluster-standard-error options as the CPB – so
compare_models() can place the CPB next to its baselines on
one footing (identical degrees of freedom, log-likelihood, and
proper-score accounting).
cpb_fit <- cpb(y ~ x, data = d, truncated = FALSE, se = "none")
compare_models(
CPB = cpb_fit,
Poisson = count_reg(y ~ x, data = d, family = "poisson"),
NB = count_reg(y ~ x, data = d, family = "negbin"),
`COM-Poisson`= count_reg(y ~ x, data = d, family = "compois")
)
#> df logLik AIC BIC logscore rps
#> CPB 3 -769.2026 1544.405 1556.380 1.923006 0.9642624
#> COM-Poisson 3 -772.0096 1550.019 1561.994 1.930024 0.9631112
#> Poisson 2 -798.9009 1601.802 1609.785 1.997252 0.9855061
#> NB 3 -798.9011 1603.802 1615.777 1.997253 0.9855063On underdispersed data the negative binomial collapses onto the
Poisson, while the CPB and COM-Poisson capture the compression and win
on AIC and the proper scores. The correlated-random-effects device
(mundlak()) and matching
d/p/q/r functions
(e.g. rcompois(), dcpb()) round out the
family.
Many count outcomes mix a participation process (most units at zero) with a tight positive count. The bundled peacekeeping panel – the number of UN operations each state contributes troops to per year – shows the package’s central move: marginally the count looks overdispersed, but conditioning on country fixed effects and benchmarking the positive counts against a zero-truncated Poisson, the at-risk process is underdispersed.
data(peacekeeping)
ud_screen(contributions ~ democracy + lgdppc + lpop + milper + factor(iso3),
data = peacekeeping, run_cpb = FALSE, run_gp = FALSE)
#>
#> === Underdispersion screen ===
#> Formula: contributions ~ democracy + lgdppc + lpop + milper + factor(iso3)
#> N=4448 mean=3.191 max=18 %zero=41.5 (unconditional var/mean=4.56)
#>
#> MARGINAL verdict: OVERDISPERSED
#> Pearson=1.086 prop.slope=0.191 (p=<2e-16)
#> NB vs Poisson LR = 18.1 (p= 1e-05 ; sig => overdispersion)
#>
#> AT-RISK (y>0) verdict:UNDERDISPERSED [n_pos=2602]
#> ZTP-Pearson = 0.841 (underdispersed if < 0.955, the calibrated 5% threshold)
#>
#> Model comparison (log-lik):
#> Poisson NB GP COMP CPB
#> -6645.816 -6636.765 NA NA NAThat is the case for a two-part model with an underdispersed
intensity. hurdle_cpb() joins a participation logit to a
zero-truncated CPB, and zi_cpb() fits the structural-zero
mixture; zi_test() and compare_models()
adjudicate between them.
z <- rnorm(n)
yh <- rhurdle_cpb(n, lambda = exp(1.2 + 0.3 * x), alpha = 0.5,
p = plogis(0.3 + 0.8 * z))
dh <- data.frame(y = yh, x = x, z = z)
h <- hurdle_cpb(y ~ x, data = dh, participation = ~ z)
zi <- zi_cpb(y ~ x, data = dh, zero = ~ z)
compare_models(hurdle = h, mixture = zi)
#> df logLik AIC BIC logscore rps
#> hurdle 5 -609.9920 1229.984 1249.941 1.524980 1.039177
#> mixture 5 -610.1172 1230.234 1250.192 1.525293 1.040127The hurdle’s first_difference() separates the extensive
and intensive margins exactly – which channel a covariate moves, not
just the blended marginal effect. One practical note: because the hurdle
factorizes, its participation stage is an ordinary logistic regression;
if a dummy-heavy participation equation separates, fit that stage with a
dedicated bias-reduction package (logistf,
brglm2) alongside this package’s zero-truncated
intensity.
The concentrated fixed-effects dispersion estimate carries the
incidental- parameters bias of order 1/T: with few observations per
unit, alpha is biased downward (the panel looks
more underdispersed than it is). bias_correct = "jackknife"
removes the leading bias term by the split-panel jackknife of Dhaene and
Jochmans (2015), refitting on each unit’s temporal halves.
short <- do.call(rbind, lapply(1:30, function(i) {
xx <- rnorm(8); NN <- pmax(round(exp(1.0 + rnorm(1, 0, 0.4) + 0.3 * xx) / 0.5), 1)
data.frame(unit = i, x = xx, y = rbinom(8, NN, 0.5))
}))
ml <- cpb_fe(y ~ x, data = short, fe = "unit")
jk <- cpb_fe(y ~ x, data = short, fe = "unit", bias_correct = "jackknife")
c(ml = ml$alpha, jackknife = jk$alpha) # truth is 0.5; ML is biased downward
#> ml jackknife
#> 0.4348484 0.5288745The correction is only valid when the two half-panels estimate the same parameter (the method’s time-homogeneity requirement), so it carries a validity gate: the panel is also split cross-sectionally by units – a placebo that is exchangeable under any time pattern – and if the temporal halves disagree beyond that placebo noise, the correction is refused with a warning naming the failed assumption and the maximum-likelihood fit is returned. On a trending or regime-changing panel, the refusal is the correct answer. The gate is deliberately powered over sized: in calibration it refuses about 9% of genuinely homogeneous panels (you keep the ordinary ML fit) while catching 98% of dispersion regime changes and all smooth unmodeled trends.
Every fitted model in the package has a simulate()
method, so the whole family plugs into DHARMa’s
simulated-residual diagnostics.
sims <- simulate(h, nsim = 100, seed = 1)
res <- DHARMa::createDHARMa(simulatedResponse = as.matrix(sims),
observedResponse = dh$y,
fittedPredictedResponse = fitted(h),
integerResponse = TRUE)
plot(res)