The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
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)
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.