---
title: "Detecting and Modeling Underdispersed Counts"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Detecting and Modeling Underdispersed Counts}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4)
set.seed(1)
```

`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.

```{r setup}
library(underdisp)
```

## Simulating an underdispersed count

We generate a count with a conditional variance-to-mean ratio of about one half.

```{r sim}
n <- 400
x <- rnorm(n)
N <- pmax(round(exp(1.6 + 0.5 * x) / 0.5), 1)
y <- rbinom(n, N, 0.5)
d <- data.frame(y = y, x = x)
c(mean = mean(y), var = var(y), ratio = var(y) / mean(y))
```

## Screening

`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).

```{r screen}
ud_screen(y ~ x, data = d, run_cpb = FALSE)
```

## Fitting the continuous parameter binomial

`cpb()` fits the CPB, with `truncated = TRUE` for the common case in which
underdispersion lives among the positive counts of a zero-inflated outcome.

```{r fit}
fit <- cpb(y ~ x, data = d[d$y > 0, ], se = "none")
summary(fit)
```

The dispersion parameter `alpha` summarizes the compression, and each observation
carries an implied ceiling `lambda / (1 - alpha)`. The CPB is a binomial with a
non-integer number of trials, renormalized over its finite support, so its mean
is not exactly the rate `exp(x'b)`: `fitted()` and
`predict(type = "response")` report the exact mean of the fitted distribution
(the conditional mean given `Y >= 1` for a zero-truncated fit), and
`predict(type = "rate")` the rate.

## Quantities of interest

Predicted means, the implied ceiling, and first differences are available for
user-specified covariate profiles.

```{r qoi}
predict(fit, newdata = data.frame(x = c(-1, 0, 1)), type = "response")
implied_ceiling(fit, newdata = data.frame(x = 0))
first_difference(fit, "x", from = -1, to = 1)
```

## Testing equidispersion

The Poisson is the equidispersed member of every family in the package (`alpha
= 1` for the CPB, `delta = 1` for the GEC, `nu = 1` for the COM-Poisson, and so
on), so `dispersion_test()` tests a fitted model's dispersion parameter against
its Poisson value by a likelihood-ratio test on the fitted design, with
fixed effects carried into the null. Where the Poisson value sits on the
boundary of the parameter space (the CPB, the negative binomial) the test is
one-sided, and its p-value comes by default from a parametric bootstrap under
the fitted Poisson, because the asymptotic Self--Liang mixture over-rejects in
finite samples; `B = 0` gives the quick asymptotic value shown here. Elsewhere
the alternative can be two-sided or directional. For a plain Poisson fit the Cameron--Trivedi
auxiliary regression is available as `method = "auxiliary"`.

```{r dtest}
dispersion_test(fit, B = 0)   # asymptotic value; the default is the bootstrap p-value
dispersion_test(count_reg(y ~ x, data = d, family = "poisson"),
                method = "auxiliary", alternative = "under")
```

## Bootstrap inference

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. That interval is
first-order by default: the estimate of `alpha` is biased toward zero (the
log-ceiling coefficients behave like endpoint parameters), so it covers about
0.90 to 0.94 in simulations, with its misses on the upper side.
`calibrate_alpha()` is the opt-in remedy: it stores the parametric-bootstrap
distribution of the signed root of the profile likelihood ratio in the fit
(`fit <- calibrate_alpha(fit, B = 199, cores = 4, seed = 1)`, about a third of
a fit per replicate), after which `alpha_confint()`, `confint()`,
`implied_ceiling()` and `summary()` report the calibrated interval. Both
intervals are model-based and assume independent observations. `cores =` runs
the replicates on a socket cluster seeded from the session.

```{r boot}
fit_b <- cpb(y ~ x, data = d[d$y > 0, ], se = "bootstrap", B = 49)   # a small B keeps the vignette quick
summary(fit_b)
irr(fit_b)             # rate ratios with percentile intervals
confint(fit_b)         # coefficients (percentile) and alpha (first-order profile likelihood)
```

## The free-dispersion GEC

`gec()` fits King's generalized event count (Katz) model, whose dispersion
`delta` 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.

```{r gec}
gec(y ~ x, data = d, se = "none")                                  # delta ~ 0.5
gec(y ~ x, data = data.frame(y = rpois(n, exp(1 + 0.4 * x)), x = x),
    se = "none")                                                   # delta ~ 1
```

The GEC carries the same zero-truncated, hurdle (`hurdle_gec()`), zero-inflated
(`zi_gec()`), and fixed-effects (`gec_fe()`) variants as the CPB.

## High-dimensional fixed effects

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.

```{r fe}
panel <- do.call(rbind, lapply(1:30, 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))
}))
fe_fit <- cpb_fe(y ~ x, data = panel, fe = "unit")
fe_fit
dispersion_test(fe_fit, B = 0)   # the statistic against a Poisson with the same unit effects
```

With unit fixed effects the asymptotic distribution of that statistic does not
apply, so its p-value comes from a parametric bootstrap under the fitted
Poisson. Every replicate refits the model, which makes the bootstrap about 199
times as slow as the fit; it is not run here:

```{r fe-boot, eval = FALSE}
dispersion_test(fe_fit, cores = 2)   # parametric-bootstrap p-value, 199 replicates
```

## Comparing the family

`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.

```{r family}
compare_dispersion(y ~ x, data = d)$table
```

## The wider underdispersed family

Underdispersion has more than one mechanism, and `count_reg()` fits the
distributions that express them through one interface: the classical
COM-Poisson (`"compois"`, rate-parameterized) and Huang's mean-parameterized
COM-Poisson (`"mpcmp"`, a soft tail), the Consul--Jain generalized Poisson
(`"genpois"`, constant variance-to-mean ratio, finite support when
underdispersed), Winkelmann's gamma-count (`"gammacount"`, events with more
regular timing than a Poisson process), and Efron's double Poisson
(`"doublepois"`). Each inherits fixed effects, zero-truncation, hurdle and
zero-inflated forms, offsets, frequency weights, robust and cluster-robust
standard errors, and every method in the package, so `compare_models()` places
them next to the CPB on one footing.

```{r wider}
dw <- d[1:150, ]                         # part of the sample keeps the seven fits quick
fits <- list(
  CPB          = cpb(y ~ x, data = dw, truncated = FALSE, se = "none"),
  Poisson      = count_reg(y ~ x, data = dw, family = "poisson"),
  NB           = count_reg(y ~ x, data = dw, family = "negbin"),
  `COM-Poisson`= count_reg(y ~ x, data = dw, family = "mpcmp"),
  GenPoisson   = count_reg(y ~ x, data = dw, family = "genpois"),
  GammaCount   = count_reg(y ~ x, data = dw, family = "gammacount"),
  DoublePois   = count_reg(y ~ x, data = dw, family = "doublepois")
)
do.call(compare_models, fits)
```

On underdispersed data the negative binomial collapses onto the Poisson, while
the underdispersed families capture the compression and win on AIC and the
proper scores. The families differ in how their dispersion moves with the mean,
and `dispersion_profile()` sets the empirical conditional variance-to-mean
ratio, by bins of the fitted mean, against the curve each family implies:

```{r profile}
dispersion_profile(CPB = fits$CPB, NB = fits$NB, GammaCount = fits$GammaCount,
                   GenPoisson = fits$GenPoisson)
```

The correlated-random-effects device (`mundlak()`) and matching
`d`/`p`/`q`/`r` functions (`dcpb()`, `dgec()`, `dcompois()`, `dgammacount()`,
`dgenpois()`, `ddoublepois()`) round out the family.

## Excess zeros: hurdle and zero-inflated models

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.

```{r pkscreen}
data(peacekeeping)
ud_screen(contributions ~ democracy + lgdppc + lpop + milper + factor(iso3),
          data = peacekeeping, run_cpb = FALSE, run_gp = FALSE)
```

That 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.

```{r hurdle}
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)
```

The 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.

## Weights

Every estimator takes frequency `weights` (a weight of `w` is equivalent to
`w` copies of the row, in the likelihood, the information, and the bootstrap),
so tabulated data fit exactly as the expanded data would:

```{r weights}
w <- sample(1:3, n, TRUE)
c(weighted = count_reg(y ~ x, data = d, family = "gammacount", weights = w)$loglik,
  expanded = count_reg(y ~ x, data = d[rep(seq_len(n), w), ], family = "gammacount")$loglik)
```

## Short panels: bias-corrected fixed effects

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.

```{r jackknife}
short <- do.call(rbind, lapply(1:16, 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
```

The 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 7% of
genuinely homogeneous panels (you keep the ordinary ML fit) while catching 98%
of dispersion regime changes and all smooth unmodeled trends.

## Simulated-residual diagnostics with DHARMa

Every fitted model in the package has a `simulate()` method, so the whole
family plugs into `DHARMa`'s simulated-residual diagnostics.

```{r dharma, eval = requireNamespace("DHARMa", quietly = TRUE)}
sims <- simulate(h, nsim = 100, seed = 1)
res <- DHARMa::createDHARMa(simulatedResponse = as.matrix(sims),
                            observedResponse  = dh$y,
                            fittedPredictedResponse = fitted(h),
                            integerResponse = TRUE)
plot(res)
```