Package {underdisp}

Type:	Package
Title:	Diagnostics and Models for Underdispersed Count Data
Version:	0.2.0
Description:	Tools for detecting and modeling underdispersion in count data (conditional variance below the conditional mean), the case the Poisson and negative binomial defaults cannot represent. Provides a screening diagnostic that benchmarks at-risk dispersion against a zero-truncated Poisson, regression-adjusted tests of equidispersion, and a dispersion profile that compares the variance-to-mean curves of competing families against the data; the continuous parameter binomial (CPB) and generalized event count (Katz) regressions with zero-truncated, hurdle, and zero-inflated forms and high-dimensional fixed effects with a split-panel jackknife bias correction; matched Poisson, negative binomial, COM-Poisson (rate- and mean-parameterized), generalized Poisson, gamma-count, double Poisson, hyper-Poisson, and exponentially weighted Poisson regressions through the same interface, with frequency weights, offsets, and analytic, robust, and cluster-robust standard errors; bootstrap and profile-likelihood inference; proper scoring rules, rootograms, PIT histograms, and simulation methods; and quantities of interest including predicted distributions, the implied ceiling, rate ratios, and first differences with an extensive/intensive decomposition. The likelihoods are implemented in C++.
License:	GPL-3
Encoding:	UTF-8
LazyData:	true
Depends:	R (≥ 4.0)
Imports:	Rcpp, stats, MASS, VGAM, graphics, methods, numDeriv, parallel
LinkingTo:	Rcpp
Suggests:	sandwich, pscl, DHARMa, testthat (≥ 3.0.0), knitr, rmarkdown, broom, modelsummary, texreg, gamlss.dist, rmutil, glmmTMB, DGLMExtPois, ewp
VignetteBuilder:	knitr
RoxygenNote:	7.3.3
URL:	https://CRAN.R-project.org/package=underdisp, https://github.com/bagozzib/underdisp, https://bagozzib.github.io/underdisp/
BugReports:	https://github.com/bagozzib/underdisp/issues
Config/testthat/edition:	3
NeedsCompilation:	yes
Packaged:	2026-11-04 15:54:03 UTC; bagoz
Author:	Benjamin E. Bagozzi [aut, cre]
Maintainer:	Benjamin E. Bagozzi <bagozzib@udel.edu>
Repository:	CRAN
Date/Publication:	2026-11-04 18:40:02 UTC

underdisp: Diagnostics and Models for Underdispersed Count Data

Description

Detect and model underdispersion (conditional variance below the conditional mean) in count data.

Screening and diagnostics

ud_screen() (the zero-truncated-Poisson at-risk screen with calibrated or parametric-bootstrap thresholds), dispersion_test() (regression-adjusted tests of equidispersion), dispersion_profile() (the conditional variance-to-mean curve against each family's implied curve), compare_dispersion(), zi_test(), rootogram(), pit_hist().

Estimators

The hard-ceiling family cpb() / cpb_fe() and the free-dispersion family gec() / gec_fe(), with hurdle and zero-inflated forms (hurdle_cpb(), zi_cpb(), hurdle_gec(), zi_gec()); the matched count families through one interface, count_reg(), hurdle_count(), zi_count() (Poisson, negative binomial, COM-Poisson in the rate and the mean parameterization, generalized Poisson, gamma-count, double Poisson, hyper-Poisson, exponentially weighted Poisson), all with offsets and frequency weights, and with fixed effects except in the count equation of the GEC's two-part forms (which takes dummies in its formula); the CPB and GEC classes take bootstrap inference, the count families analytic, robust or clustered errors.

Comparison and quantities of interest

compare_models(), score(), cv_score(); predict(), implied_ceiling(), ceiling_exceedance(), irr(), first_difference() (with the extensive/intensive decomposition for the two-part models), ame(), confint() for every class, simulate() for DHARMa diagnostics, and broom / texreg / modelsummary support.

Author(s)

Maintainer: Benjamin E. Bagozzi bagozzib@udel.edu (ORCID)

See Also

Useful links:

• https://CRAN.R-project.org/package=underdisp

• https://github.com/bagozzib/underdisp

• https://bagozzib.github.io/underdisp/

• Report bugs at https://github.com/bagozzib/underdisp/issues

Profile-likelihood interval for the dispersion parameter alpha

Description

The interval inverts the likelihood-ratio test of alpha, re-maximizing the coefficients at each value. By default the cut is the chi-square one (a first-order interval). On King's support it covers between 0.940 and 0.956 in simulations from the CPB. On the floor support it covers about 0.90 to 0.94, with most misses on the upper side, because the estimate of alpha is biased toward zero there (calibrate_alpha() explains the mechanism). A fit that went through calibrate_alpha() gets the interval calibrated by parametric bootstrap instead. Either interval is model-based: it assumes the CPB and independent observations.

Usage

alpha_confint(object, level = 0.95)

Arguments

Value

A length-2 numeric vector (lower, upper) with attributes alpha (the point estimate), method (first-order or calibrated) and boundary (TRUE when the profile does not fall to the cut on the lower side: on King's support, whose trace runs down to alpha = 1e-6, the lower limit is then 0; on the floor support, whose trace stops at alpha = 0.005, it is that bound).

Examples

set.seed(7); x <- rnorm(200)
N <- pmax(round(exp(1.5 + 0.4 * x) / 0.5), 1); y <- rbinom(200, N, 0.5)
fit <- cpb(y ~ x, data.frame(y = y, x = x)[y > 0, ], se = "none")
alpha_confint(fit)

Average marginal effects

Description

The average, over the rows of newdata, of the change in a predicted quantity when one variable changes: its derivative (a numeric variable, with from and to left out) or the contrast between two of its values (from and to, for a variable of any type). Every row keeps its own values of the other variables, and the variable enters through the model's formula, so interactions and transformations of it (x:z, I(x^2), log(x)) are handled. first_difference() instead moves one model-matrix column at a single reference profile.

Usage

ame(
  object,
  variable,
  newdata,
  from = NULL,
  to = NULL,
  type = "response",
  eps = NULL,
  level = 0.95,
  offset = NULL,
  weights = NULL,
  B = 0,
  cores = 1L
)

Arguments

Details

The interval follows the rule of first_difference(). With stored bootstrap draws, the average is computed at every draw, jointly over the coefficients, the dispersion parameter and a fixed-effects fit's reference unit (the draws of a fixed-effects fit recentered on its estimate), and the percentile interval is taken. Otherwise it uses the delta method over the fit's covariance, and otherwise there is none. A row at a fixed-effects unit's own effect, or at its own zero-equation dummy, depends on a parameter the draws do not carry, so an average over such rows, the sample's average partial effect, has no interval of this kind (B below gives it one). Without the unit column in newdata every row enters at the reference unit, the mean of the identified units' effects on the log scale (a typical, geometric-mean unit), whose draws the fit carries: that is the average effect for such a unit over the rows' covariates, not the sample's average partial effect, and the two can differ substantially when the unit effects are dispersed.

A fixed-effects count family (count_reg() with fe) carries its unit dummies in its covariance, so its average partial effect has a delta-method interval. The slope of a fit on the floor support (support = "floor") is refused, since its mean jumps at integer ceilings; a contrast is available.

B > 0 gives the average a refit-bootstrap interval instead, for any class: the model is refit on B resampled copies of newdata, which must then be the data the model was fit to, resampling whole clusters of the fit's cluster, else whole units of its first fixed-effects column, else rows (a unit drawn twice becomes two units). Each replicate averages over its own rows at its own unit effects, and the percentile interval is taken. Resampling the units treats them as a sample, so this is an interval for a fixed-effects fit's average partial effect in the population of units they represent, which carries the variability of which units are drawn: with few units, or a few with large effects, it is wide and asymmetric. Each refit is the fit's own call with its data replaced and its standard errors off (on King's support, every replicate of a jackknife-corrected fit is corrected without re-running the validity gate, as in cpb_fe()'s own bootstrap, and none is when the gate refused the fit's correction); the call's other arguments are evaluated where ame() is called, as update() does, and an offset or weights given to the fit as one value per row follow the resampled rows. The run takes about B times one fit (see cores).

Value

A data frame of class "ud_ame" with term, contrast ("dY/dx", or "to - from"), type, estimate, lower, upper, method (the interval's source) and n, the rows averaged (a row with a missing value is left out).

See Also

first_difference(), predict.cpb() (se.fit, interval).

Examples

set.seed(1); n <- 400; x <- rnorm(n); g <- factor(sample(c("a", "b"), n, TRUE))
d <- data.frame(y = rcpb(n, exp(1 + 0.4 * x + 0.3 * (g == "b")), 0.5), x = x, g = g)
fit <- cpb(y ~ x * g, data = d, truncated = FALSE, se = "bootstrap", B = 100)
ame(fit, "x", d)
ame(fit, "g", d, from = "a", to = "b")

## a fixed-effects fit's average partial effect over its own units
pd <- do.call(rbind, lapply(1:30, function(i) {
  x <- rnorm(8); data.frame(unit = i, x = x, y = rcpb(8, exp(rnorm(1, 0, 0.5) + 0.4 * x), 0.5))
}))
ffit <- cpb_fe(y ~ x, data = pd, fe = "unit")
ame(ffit, "x", pd, B = 100)

Augment data with CPB fitted values (broom method)

Description

Augment data with CPB fitted values (broom method)

Usage

## S3 method for class 'cpb'
augment(x, data = NULL, ...)

Arguments

Value

A data frame with fitted values, response residuals, and the implied per-observation ceiling (appended to data when it is given).

Calibrate the interval for the CPB dispersion parameter

Description

An opt-in step that stores, in the fit, the parametric-bootstrap distribution of the signed root of the profile likelihood ratio for alpha. Once it is stored, alpha_confint(), confint.cpb(), implied_ceiling() and summary() report the calibrated interval instead of the first-order one. It is the remedy for floor-support fits (support = "floor"). On King's support, the default, the first-order interval covers close to its level, and the calibration is an optional refinement.

Usage

calibrate_alpha(object, B = 199, cores = 1L, seed = NULL, ...)

## S3 method for class 'cpb'
calibrate_alpha(object, B = 199, cores = 1L, seed = NULL, ...)

## S3 method for class 'hurdle_cpb'
calibrate_alpha(object, B = 199, cores = 1L, seed = NULL, ...)

## Default S3 method:
calibrate_alpha(object, B = 199, cores = 1L, seed = NULL, ...)

Arguments

Details

The first-order interval cuts the profile log-likelihood qchisq(level, 1) / 2 below its maximum on both sides. On King's support that cut is close to calibrated, because the probabilities change continuously with the parameters. The simulations from the CPB cover nine designs (n from 200 to 1,000, true alpha 0.3, 0.5 and 0.7). There the first-order 95% interval covers between 0.940 and 0.956, 0.947 on average, and more of its misses fall above the upper limit than below the lower one. With B = 39 replicates, the calibrated interval covers 0.948 on average in the same simulations. It is about 12% longer, and its misses are more nearly balanced.

On the floor support the cut is not calibrated. The support ⁠0, ..., floor(lambda / (1 - alpha))⁠ moves with the parameters, the log-ceiling coefficients behave like endpoint parameters, and alpha-hat inherits a bias toward zero from them. Across the same nine designs the first-order 95% interval covers between 0.90 and 0.94, and most of its misses have the true alpha above the upper limit.

The calibration simulates B responses from the fitted model (the fitted rates, alpha-hat, the same truncation and offset), refits each one, and records the signed root ⁠r* = sign(alpha* - alpha-hat) sqrt(2 (l* - l*_p(alpha-hat)))⁠. The interval is then ⁠{alpha : q_lo <= r(alpha) <= q_hi}⁠, with q_lo and q_hi the order statistics of the ⁠r*⁠ at ranks (B + 1)(1 - level)/2 from each end (the 5th smallest and 5th largest of 199 at the 95% level), so the upper limit sits q_lo^2 / 2 and the lower limit q_hi^2 / 2 below the maximum of the profile. On the floor support, in the same simulations, the calibrated 95% interval covers 0.92 to 0.96, with its misses balanced between the two sides. On the floor support, the first-order interval falls short only when the design has many distinct covariate patterns (the endpoint effect needs many distinct ceilings): with an intercept only, or a covariate taking ten or fewer values, it covers at or above its level, and the calibration, which stays close to the nominal level there too (0.95 with a binary covariate), is not needed.

The responses are all drawn first, on the calling process, and the refits use no random numbers, so the result depends on seed but not on cores. With B = 199 the two cuts carry simulation error (a standard deviation of about 0.2 in q_lo and q_hi), which moves a limit by a fraction of its distance from the estimate; B = 999 reduces it and is needed for a 99% interval. Each replicate costs about a third of the original fit.

The interval is model-based. It assumes the CPB and independent observations, and it is not cluster-robust: when the fit was given cluster, a message says so. Under misspecification alpha has no fixed target (on data with a soft upper tail its estimate rises with the sample size), and no interval computed under the fitted CPB repairs that.

Calibration is refused, with the reason, where it has nothing to work with: a fit whose ceiling reaches max.support (the guard, not the data, bounds alpha), weights that are not whole numbers, and a fit saved by version 0.1.0. Whole-number weights are treated as frequency weights: a row of weight w contributes w independent simulated responses.

Value

object, with the component alpha.cal: a list holding the signed roots r (NA for a replicate that failed to fit), B, B_ok and seed.

See Also

alpha_confint(), dispersion_test()

Examples

set.seed(7); x <- rnorm(100)
d <- data.frame(y = rcpb(100, exp(1.5 + 0.4 * x), 0.5, truncated = TRUE), x = x)
fit <- cpb(y ~ x, d, se = "none")
alpha_confint(fit)                       # first-order
## 19 replicates keep the example short and can serve a 90% interval; use 199 or more
fit <- calibrate_alpha(fit, B = 19, seed = 1)
alpha_confint(fit, level = 0.90)         # calibrated

Held-out ceiling exceedance of a CPB fit

Description

Whether counts respect the CPB's ceiling can be checked, but not in sample: its likelihood is minus infinity for any count above the bound, so every feasible fit respects the bound on the data it was fit to. The check is one-sided: it detects a ceiling that is too low, where held-out counts reach it, and has no power against one that is too high (a fitted ceiling above a true bound is never exceeded), nor where counts stay well below the ceiling (see implied_ceiling()). ceiling_exceedance() takes held-out observations and reports how many of their counts exceed the largest count the fitted model allows at their covariates (max_count of implied_ceiling()); for a two-part model the bound is the count component's, and a fixed-effects unit seen in training enters at its own effect (an unseen one at the mean effect), as in score(). A breached ceiling also zeroes the held-out likelihood, which score() records with its 1e-12 floor.

Usage

ceiling_exceedance(fit, newdata, weights = NULL, offset = NULL)

Arguments

Value

A one-row data frame with n (held-out rows, weighted), exceed (rows above the bound, weighted), and share.

See Also

implied_ceiling(), score()

Examples

set.seed(7); x <- rnorm(400); y <- rcpb(400, exp(1.4 + 0.3 * x), 0.5)
d <- data.frame(y = y, x = x)
fit <- cpb(y ~ x, d[1:300, ], truncated = FALSE, se = "none")
ceiling_exceedance(fit, d[301:400, ])

Compare count-model dispersion across the model family

Description

Fits the Poisson and negative binomial defaults alongside the underdispersed models—the soft-tail Conway–Maxwell–Poisson (fit natively, no external dependency), the hard-ceiling CPB, and the GEC, which spans both directions—and returns a comparison on both fit (log-likelihood, AIC, BIC) and calibration (mean logarithmic score and ranked probability score, lower is better, plus the fitted share of zeros against the observed share). Whether the tail is better described by an accelerated decay (COM-Poisson) or a hard ceiling (CPB) is a question for the data's fit and calibration, not an assumption; the scoring rules are the calibration counterpart to the information criteria. Either family's better fit describes the distribution, not the process behind it. The CPB's bound can be checked only on held-out counts, with ceiling_exceedance(), and only against a ceiling that is too low: on the data it was fit to every feasible CPB respects its bound, since its likelihood forbids exceedance.

Usage

compare_dispersion(
  formula,
  data,
  max.support = NULL,
  hurdle = FALSE,
  zi = FALSE,
  support = c("ceiling", "floor")
)

Arguments

Details

With hurdle = TRUE and zeros present, a hurdle-CPB is added, so a zero-inflated underdispersed process can be compared to the single-equation models on the same footing.

Value

A list with table (a data frame of df, logLik, AIC, BIC, logscore, rps, zero_fit, and converged per model), obs_zero (the observed zero share), nu (the COM-Poisson dispersion, ⁠>1⁠ = underdispersion, NA if the COM-Poisson fit failed), alpha (the CPB shape parameter), gec_delta (the GEC's dispersion, NA if its fit failed), and cpb_ok (FALSE if the single-equation CPB could not satisfy its feasibility constraint on the data, e.g. under heavy zero-inflation with a wide count range: the CPB's support cannot hold the counts, which another family, or a hurdle or zero-inflated model, may fit).

See Also

ud_screen(), cpb(), score()

Examples

set.seed(1); x <- rnorm(120)
N <- pmax(round(exp(1.6 + 0.4 * x) / 0.5), 1); y <- rbinom(120, N, 0.5)
compare_dispersion(y ~ x, data = data.frame(y = y, x = x))$table

Compare fitted underdispersed-count models

Description

Compares fitted models from this package — cpb, cpb_fe, hurdle_cpb, zi_cpb, gec, gec_fe, hurdle_gec, zi_gec, count_reg, hurdle_count, and zi_count, in any combination of fixed effects and robust/clustered standard errors — on information criteria and, where a predicted distribution is available, on proper scores. Unlike zi_test() this is not a hypothesis test: the models need not be nested, so it is the appropriate tool for the non-nested hurdle-versus-mixture choice. Passing an object that is not a fitted model from this package is an error, and models fit on different data raise a warning.

Usage

compare_models(...)

Arguments

Value

A data frame with df, logLik, AIC, BIC, and (where the predicted distribution is available) logscore and rps, one row per model, ordered by AIC, and the flags truncated (a zero-truncated likelihood, conditional on Y >= 1) and separated (a zero-inflation or participation equation the data do not bound, whose log-likelihood is a supremum along a divergent ridge; see zi_cpb(), Separation). The function warns when the models mix truncated and untruncated likelihoods, whose log-likelihoods are of different events, or when a model is separated.

Degrees of freedom

Each model's df is its own (logLik()'s). The hurdle and zero-inflated fits do not count the parameters of a unit with no positive count: with unit fixed effects in the participation or zero equation (part_fe, zero_fe) or in the count equation (fe; dummies written into a formula are covariates to the fit and are counted), such a unit's dummy and its count effect sit at their limits, where every one of its rows has probability one, so it adds nothing to the likelihood. As fixest leaves such units out of a fixed-effects fit, these parameters are left out of df, so that a hurdle and a zero-inflated model fit to the same panel are charged alike. The fits record the units (part_zero_units in the hurdles; zero_units and zero_fe_zero_units in zi_cpb(); fe_zero_units and zero_fe_zero_units in zi_count(); zero_fe_zero_units in zi_gec()), and their print methods give the number. In new data (predict(), score(), cv_score()) the hurdles give such a unit's rows the participation probability at that limit, zero, as the zero-inflated fits give them a count effect at its limit. A unit whose every count is positive has its hurdle participation dummy at the other limit, plus infinity; that dummy is counted, as the zero-inflated fits count their zero-equation dummy for the unit, and in new data the hurdles give the unit's rows the participation probability one (the fits record these units in part_one_units). The zero-inflated fits give a zero_fe unit's rows the structural-zero probability at its limit in the same way: one for a unit with no positive count, zero for a unit whose every count is positive (zero_fe_one_units). CPB fits on the floor support (support = "floor") count these parameters, and in new data evaluate these units at the coefficients the fit stopped at.

The rule is the two-part fits'. The single-equation fixed-effects fits (cpb_fe(), gec_fe(), count_reg(fe = )) count every unit's effect, a unit with no positive count included; set beside a two-part fit to the same panel, such a fit is charged two AIC points per unit with no positive count that the two-part fit is not, and compare_models() says so in a message.

See Also

zi_test(), compare_dispersion()

Examples

set.seed(1); n <- 400; x <- rnorm(n); z <- rnorm(n)
y <- rhurdle_cpb(n, exp(1.2 + 0.5 * x), 0.5, plogis(-0.2 + 0.8 * z))
d <- data.frame(y = y, x = x, z = z)
compare_models(hurdle = hurdle_cpb(y ~ x, data = d, participation = ~ z),
               zi = zi_cpb(y ~ x, data = d, zero = ~ z))

Conway–Maxwell–Poisson distribution functions

Description

Density, distribution, quantile, and random generation for the COM-Poisson with rate lambda (or mean mu) and dispersion nu (nu > 1 underdispersed, nu = 1 Poisson, nu < 1 overdispersed). Complements the estimators count_reg(..., family = "compois") (rate parameterization) and count_reg(..., family = "mpcmp") (mean parameterization).

Usage

dcompois(x, lambda, nu, log = FALSE, mu = NULL)

pcompois(q, lambda, nu, lower.tail = TRUE, log.p = FALSE, mu = NULL)

qcompois(p, lambda, nu, lower.tail = TRUE, log.p = FALSE, mu = NULL)

rcompois(n, lambda, nu, mu = NULL)

Arguments

Value

dcompois a density, pcompois a CDF, qcompois a quantile, rcompois a numeric vector of count draws. A rate whose normalizing sum the 100,000-term cap would truncate (a mean of order 10^5 or more) gives NA.

Examples

dcompois(0:5, lambda = 3, nu = 1.5)
mean(rcompois(1000, lambda = 3, nu = 1.5))

Confidence intervals for a CPB fit

Description

Coefficient intervals use the cold-multistart bootstrap percentile method. In simulations from King's CPB (n = 300) the 95% interval for a slope covers 0.91 at alpha 0.3, 0.95 at 0.5 and 0.94 at 0.7, and on the floor support 0.95. The interval for alpha is a profile-likelihood interval. It is first-order by default: on King's support it covers close to its level in simulations from the CPB, and on the floor support about 0.90 to 0.94, with most misses on the upper side (see calibrate_alpha()). It is calibrated by parametric bootstrap when the fit went through calibrate_alpha().

Usage

## S3 method for class 'cpb'
confint(object, parm, level = 0.95, ...)

Arguments

Value

A matrix of lower/upper bounds.

Confidence intervals for every model class

Description

confint() methods for the whole family, one interval rule per inference type: bootstrap percentile intervals for cpb() (coefficients; the dispersion parameter gets its profile-likelihood interval, first-order or calibrated by calibrate_alpha()) and gec() (coefficients and delta); normal-approximation intervals from the bootstrap standard errors for the fixed-effects fits (cpb_fe(), gec_fe()) and the zero-inflated mixtures (zi_cpb(), zi_gec()); Wald intervals from the analytic covariance for count_reg() and zi_count() (the dispersion parameter's interval is mapped from its estimation scale to the natural scale); and, for the hurdles, the participation model's Wald intervals (prefixed ⁠participation:⁠) stacked over the intensity model's intervals (prefixed ⁠intensity:⁠); the zero-inflated classes prefix ⁠count:⁠ and ⁠zero:⁠. The same names label vcov(), tidy(), and the texreg tables. A fit without inference errors informatively.

Usage

## S3 method for class 'gec'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'cpb_fe'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'count_reg'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'zi_count'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'hurdle_cpb'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'hurdle_gec'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'hurdle_count'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'zi_cpb'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'zi_gec'
confint(object, parm, level = 0.95, ...)

Arguments

Value

A two-column matrix of lower and upper bounds.

See Also

confint.cpb()

Matched count regressions: Poisson, negative binomial, COM-Poisson, generalized Poisson, gamma-count, double Poisson, hyper-Poisson, exponentially weighted Poisson

Description

Fits a count regression with a log link from any of the package's count families, optionally zero-truncated and/or with unit fixed effects, returning an object that compare_models(), score(), dispersion_test(), dispersion_profile(), and the broom methods treat on the same footing as a cpb() or gec() fit. The families share one interface, so a gamma-count and a COM-Poisson can be compared to a CPB with identical degrees-of-freedom, log-likelihood, proper-score, weight, and robust-standard-error accounting rather than reconciled across packages.

Usage

count_reg(
  formula,
  data,
  family = c("poisson", "negbin", "compois", "mpcmp", "genpois", "gammacount",
    "doublepois", "hyperpois", "ewp3", "ewp2"),
  truncated = FALSE,
  fe = NULL,
  offset = NULL,
  weights = NULL,
  se = c("analytic", "robust", "cluster", "none"),
  cluster = NULL
)

Arguments

Value

An object of class "count_reg". The component ⁠$theta⁠ holds the dispersion parameter on its natural scale (the negative-binomial size, the COM-Poisson \nu, the generalized-Poisson \lambda, the gamma-count \alpha, the double-Poisson \theta, the hyper-Poisson \gamma, the EWP's (\beta_1, \beta_2) or \beta; NA for the Poisson); ⁠$fitted.values⁠ and ⁠$residuals⁠ are on the mean scale (the conditional mean E(Y | Y >= 1) for a zero-truncated fit), as fitted() and residuals() return them, and ⁠$mu⁠ is the family's natural parameter ⁠exp(offset + x'b)⁠.

Families

• "poisson": the equidispersed baseline.

• "negbin": negative binomial (overdispersion only); ⁠$theta⁠ is the size.

• "compois": the classical Conway–Maxwell–Poisson in its rate parameterization, \log\lambda = x'\beta. The rate \lambda is not the mean (though predict(type = "response") and fitted() return the mean); ⁠$theta⁠ is the dispersion \nu (\nu > 1 underdispersed, \nu = 1 Poisson, \nu < 1 overdispersed).

• "mpcmp": the COM-Poisson in Huang's (2017) mean parameterization, \log\mathrm{E}(Y) = x'\beta, with the same dispersion \nu; the coefficients are effects on the log mean and rate ratios are exact.

• "genpois": the Consul–Jain generalized Poisson with constant dispersion \lambda \in (-1, 1) and \mu = \exp(x'\beta), the mean for \lambda \ge 0, with \mathrm{Var}/\mathrm{Mean} = 1/(1-\lambda)^2; \lambda < 0 is underdispersion, on a finite support that the pmf is renormalized over, and the mean then departs from \mu where that support is short (small \mu at strong underdispersion: about 1.19\mu at \lambda = -0.6, \mu = 0.5); the exact mean is reported (see genpois-distribution).

• "gammacount": Winkelmann's (1995) gamma-count renewal-process model, underdispersed when the waiting times between events are more regular than exponential (\alpha > 1); \exp(x'\beta) is the long-run event rate and the exact mean is reported (see gammacount-distribution).

• "doublepois": Efron's (1986) double Poisson with the exact normalizing constant, \theta > 1 underdispersed (see doublepois-distribution).

• "hyperpois": the hyper-Poisson of Bardwell and Crow (1964) in the mean parameterization of Sáez-Castillo and Conde-Sánchez (2013), mean \mu = \exp(x'\beta), \gamma < 1 underdispersed, with one \gamma for all observations: their constant-dispersion model hP(\mu_i, \gamma). Their model with covariates in the dispersion, hP(\mu_i, \gamma_i), is fit by DGLMExtPois::glm.hP(). Its underdispersion is bounded: as \gamma \to 0 the variance falls only to \mu - 1 for \mu \ge 1 and \mu(1 - \mu) below, so Var/Mean is at least 0.5 at a mean of 2 and 0.9 at 10 (see hyperpois-distribution).

• "ewp3" (alias "ewp"): the three-parameter exponentially weighted Poisson of Ridout and Besbeas (2004), a Poisson with centrality \lambda = \exp(x'\beta) reweighted by a two-sided exponential in the distance from it; ⁠$theta⁠ holds the two dispersion parameters (\beta_1, \beta_2), one per tail, positive values underdispersed (see ewp-distribution). "ewp2" is its two-parameter form, \beta_1 = \beta_2 = \beta. The coefficients are effects on the log centrality, not the log mean; the exact mean is reported.

The Poisson is the equidispersed member of every family, so dispersion_test() tests each family's dispersion parameter against it, and dispersion_profile() compares the families' implied variance-to-mean curves against the data.

Runtime: the Poisson, negative binomial, gamma-count, and double Poisson and generalized Poisson families fit in under a second at a few hundred rows; the COM-Poisson families, the hyper-Poisson and the EWP evaluate a normalizing sum per observation (and the mean parameterizations solve a root per observation), so they take seconds at a few hundred rows and minutes at several thousand. The EWP's sum lengthens with a heavy upper tail (a negative upper weight), and its standard errors sum over each observation's support. Rows that repeat (the same count, covariates and offset) are evaluated once with their total weight, so an intercept-only or few-category design costs as many evaluations as it has distinct rows. se = "robust" and "cluster" add numerical scores at the same cost per parameter.

References

King, G. and Roberts, M. E. (2015). How robust standard errors expose methodological problems they do not fix, and what to do about it. Political Analysis, 23(2), 159-179. Huang, A. (2017). Mean-parametrized Conway–Maxwell–Poisson regression models for dispersed counts. Statistical Modelling, 17(6), 359-380. Winkelmann, R. (1995). Duration dependence and dispersion in count-data models. Journal of Business & Economic Statistics, 13(4), 467-474. Efron, B. (1986). Double exponential families and their use in generalized linear regression. Journal of the American Statistical Association, 81(395), 709-721. Consul, P. C. and Famoye, F. (1992). Generalized Poisson regression model. Communications in Statistics – Theory and Methods, 21(1), 89-109. Shmueli, G., Minka, T. P., Kadane, J. B., Borle, S. and Boatwright, P. (2005). A useful distribution for fitting discrete data: revival of the Conway–Maxwell–Poisson distribution. Journal of the Royal Statistical Society C, 54(1), 127-142. Sellers, K. F. and Shmueli, G. (2010). A flexible regression model for count data. Annals of Applied Statistics, 4(2), 943-961. Bardwell, G. E. and Crow, E. L. (1964). A two-parameter family of hyper-Poisson distributions. Journal of the American Statistical Association, 59(305), 133-141. Sáez-Castillo, A. J. and Conde-Sánchez, A. (2013). A hyper-Poisson regression model for overdispersed and underdispersed count data. Computational Statistics & Data Analysis, 61, 148-157. Ridout, M. S. and Besbeas, P. (2004). An empirical model for underdispersed count data. Statistical Modelling, 4(1), 77-89. Gourieroux, C., Monfort, A. and Trognon, A. (1984). Pseudo maximum likelihood methods: applications to Poisson models. Econometrica, 52(3), 701-720. Wooldridge, J. M. (1999). Distribution-free estimation of some nonlinear panel data models. Journal of Econometrics, 90(1), 77-97.

See Also

cpb(), gec(), compare_models(), dispersion_test(), dispersion_profile(), hurdle_count(), zi_count()

Examples

set.seed(1); n <- 400; x <- rnorm(n)
y <- rgammacount(n, exp(1 + 0.5 * x), alpha = 2)     # regular event timing
d <- data.frame(y = y, x = x)
m <- count_reg(y ~ x, data = d, family = "gammacount")
m
compare_models(poisson = count_reg(y ~ x, d, family = "poisson"), gammacount = m)

Fit a Continuous Parameter Binomial (CPB) regression

Description

Fits the underdispersed continuous parameter binomial model of King (1989): a binomial count with a continuous number of trials, \lambda_i/(1-\alpha), and success probability 1-\alpha, so that the conditional variance is close to a fraction of the conditional mean, \mathrm{Var}(Y\mid x) \approx \alpha\,\mathrm{E}(Y\mid x) with 0 < \alpha < 1, and each observation has an endogenous ceiling, \lceil \lambda_i/(1-\alpha) \rceil. A zero-truncated variant (the default) conditions on Y \ge 1, appropriate when underdispersion lives among the positive counts of an otherwise zero-inflated outcome.

Usage

cpb(
  formula,
  data,
  truncated = TRUE,
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  offset = NULL,
  weights = NULL,
  cores = 1L,
  alpha.start = 0.5,
  max.support = NULL,
  maxit = 20000,
  reltol = 1e-08,
  support = c("ceiling", "floor")
)

Arguments

Details

The rate is modeled log-linearly, \lambda_i = \exp(x_i'\beta). The probabilities follow the Katz recursion P(Y = k+1)/P(Y = k) = (\lambda_i - k(1-\alpha))/(\alpha(k+1)) on the support 0, \ldots, \lceil \lambda_i/(1-\alpha) \rceil, King's normalization as Winkelmann, Signorino and King (1995) give it. The CPB is therefore the generalized event count model of gec() at \delta = \alpha, and the two functions return the same fit wherever the GEC's estimate of \delta is below one, given the same truncated (cpb() defaults to the zero-truncated model, gec() to the untruncated one). Where \lambda_i/(1-\alpha) is an integer the distribution is the binomial with that many trials, with mean \lambda_i and variance \alpha\lambda_i; elsewhere its mean is at least \lambda_i. In King's notation alpha is \sigma^2, the variance-to-mean ratio at an integer ceiling; Winkelmann, Signorino and King write the ceiling \lambda_i/(1-\alpha) as \alpha and 1-\alpha as p. Because the support moves with the parameters, the log-likelihood is continuous with kinks wherever an observation's ceiling crosses an integer, and it tends to the Poisson log-likelihood as \alpha \to 1; a fit that reaches that limit (\hat\alpha > 0.999) warns, and dispersion_test() gives it the boundary value of the test. The maximizer is found by a deterministic multistart, Nelder-Mead from the Poisson solution at alpha.start and seven further dispersions, the best polished by Nelder-Mead with feasibility-repaired restarts and a final BFGS pass, on covariates scaled to unit standard deviation (the coefficients are mapped back, so the fit does not depend on the covariates' units), and the reported maximum 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 neighboring solution) until it has dropped four log-likelihood units, and a trace point above the multistart's value restarts the fit from there. The same trace gives the profile interval of confint.cpb(). Unit intercepts belong in the fixed-effects estimator cpb_fe(), which solves each of them exactly (a pooled fit with many dummy columns can stop short of it).

The likelihood is not regular at the crossings. Just above an integer the top support point carries probability proportional to the ceiling's distance \epsilon past it, so the information an observation carries about its rate grows like 1/\epsilon (about 0.28/\epsilon for \log\lambda at \alpha = 0.5 and a ceiling of 3 + \epsilon, about 1.8/\epsilon at \alpha = 0.3; King 1989, footnote 7, notes that asymptotic normality needs non-integer ceilings). With continuous covariates the ceilings' fractional parts are spread and the expected information per observation is unbounded, so standard errors from the numerical Hessian are unreliable (at small alpha they fall well below the sampling variability of the estimates), and inference uses a cold-multistart bootstrap for the coefficients and a profile-likelihood interval for \alpha. In simulations from King's CPB at n = 300, the bootstrap's 95% interval for a slope covers 0.91 at alpha 0.3, 0.95 at 0.5 and 0.94 at 0.7 (0.92 at alpha 0.3 with n = 1,000); on the floor support it covers 0.95. The interval for \alpha is first-order by default and calibrated by parametric bootstrap after calibrate_alpha() (see alpha_confint()).

support = "floor" gives the package's definition before version 0.2.0, on 0, \ldots, \lfloor \lambda_i/(1-\alpha) \rfloor, whose mean is at most \lambda_i. Its log-likelihood is discontinuous: whenever a parameter move carries an observation's ceiling across an integer its normalizing constant jumps, by about (1-\alpha)^{k} at ceiling k. The fit then adds a scan of the breakpoints in alpha next to the estimate; two parameterizations of the same design (say, treatment and sum contrasts) can settle on different teeth of the surface, with log-likelihoods differing by the order of the jumps; and ceilings above max.support are infeasible, so alpha can be bounded by that guard rather than by the data.

Runtime: a fit with se = "none" takes about a second at 500 rows and four at 2,000 on one core, the multistart and the profile trace included (about four and fifteen on the floor support). The bootstrap replicates are maximized by the multistart without the profile trace (about 0.1 seconds each at 500 rows; on the floor support a replicate can therefore sit on a different tooth from the point estimate, by the order of the jumps); cores runs them in parallel.

Value

An object of class "cpb": a list with coefficients, alpha, bootstrap se.beta/vcov/ci.beta, loglik, loglik.null, fitted.values, ceiling (the observation-specific implied ceiling), the model frame pieces, and the original call.

References

King, G. (1989). Variance specification in event count models: From restrictive assumptions to a generalized estimator. American Journal of Political Science, 33(3), 762-784. Winkelmann, R., Signorino, C. S. and King, G. (1995). A correction for an underdispersed event count probability distribution. Political Analysis, 5, 215-228. King, G. and Roberts, M. E. (2015). How robust standard errors expose methodological problems they do not fix, and what to do about it. Political Analysis, 23(2), 159-179.

See Also

ud_screen(), confint.cpb(), predict.cpb(), implied_ceiling(), dispersion_test()

Examples

set.seed(1)
n <- 200; x <- rnorm(n)
N <- pmax(round(exp(1.6 + 0.5 * x) / 0.5), 1)
y <- rbinom(n, N, 0.5)                 # underdispersed (var/mean approx 0.5)
d <- data.frame(y = y, x = x)
fit <- cpb(y ~ x, data = d[d$y > 0, ], se = "none")
summary(fit)

Continuous parameter binomial distribution functions

Description

Density, distribution function, and quantile function for the continuous parameter binomial (CPB), and its zero-truncated form. The CPB has a hard ceiling: King's CPB (the default) puts no probability above \lceil \lambda/(1-\alpha) \rceil, the floor support (support = "floor", the definition before version 0.2.0) none above \lfloor \lambda/(1-\alpha) \rfloor. Complements the simulator rcpb().

Usage

dcpb(
  x,
  lambda,
  alpha,
  truncated = FALSE,
  log = FALSE,
  support = c("ceiling", "floor")
)

pcpb(
  q,
  lambda,
  alpha,
  truncated = FALSE,
  lower.tail = TRUE,
  log.p = FALSE,
  support = c("ceiling", "floor")
)

qcpb(
  p,
  lambda,
  alpha,
  truncated = FALSE,
  lower.tail = TRUE,
  log.p = FALSE,
  support = c("ceiling", "floor")
)

Arguments

Value

dcpb a density, pcpb a distribution function, qcpb a quantile.

See Also

rcpb()

Examples

dcpb(0:5, lambda = 3, alpha = 0.5)
pcpb(3, lambda = 3, alpha = 0.5)
qcpb(0.9, lambda = 3, alpha = 0.5)

CPB regression with high-dimensional unit fixed effects

Description

Fits the CPB with a full set of unit fixed effects by concentrating (profiling) out the unit intercepts. Each unit's intercept is solved by a one-dimensional inner maximization, so the outer optimizer handles only the covariate coefficients and the dispersion parameter. This avoids an explicit unit-dummy design matrix and scales to thousands of units, where dummy-based fitting exhausts memory. The concentrated log-likelihood equals the full-dummy log-likelihood to the tolerance of the inner search, which is solved exactly. On King's support (the default) each unit's objective is continuous in its intercept: when the intercept carries observation t's ceiling across an integer the support point that enters carries mass growing from zero, so the crossing is a kink, and the objective falls to minus infinity at the feasibility floor, where some count would exceed \lceil \lambda_{it}/(1-\alpha) \rceil. As for the GEC (the same distribution at \delta = \alpha; see gec_fe()), the smooth part is unimodal and the kinks only lower the slope, so the inner search is a golden section on a bracket around the unit's Poisson intercept, expanded while the maximum sits at an edge; a maximum at a kink is tracked in the analytic gradient of the concentrated likelihood, which the outer BFGS uses. The outer optimizer is a deterministic multistart: BFGS with that gradient from the Poisson slopes and two dispersion starts, a Nelder-Mead polish, and two perturbed restarts.

Usage

cpb_fe(
  formula,
  data,
  fe,
  truncated = FALSE,
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  offset = NULL,
  weights = NULL,
  cores = 1L,
  max.support = NULL,
  inner_it = 30L,
  maxit = 3000L,
  reltol = 1e-07,
  bias_correct = c("none", "jackknife"),
  time = NULL,
  support = c("ceiling", "floor")
)

Arguments

Details

On the floor support (support = "floor", the definition before version 0.2.0) each unit's objective is a saw-tooth in its intercept: whenever the intercept crosses \log(k(1-\alpha)) - o_{it} for an integer k the support of observation t gains the point k, its normalizing constant grows, and the objective drops by about (1-\alpha)^k; between such breakpoints it is smooth and unimodal. Its supremum is therefore attained at an interior critical point of one tooth or as the left limit at a breakpoint, and the inner search evaluates every breakpoint within a window of the envelope's peak (located on a coarse grid around the unit's Poisson intercept) and refines inside whole teeth by golden section; every evaluation of the concentrated likelihood searches each unit afresh, so the objective is a function of the parameters alone. Solving the inner problem to its supremum makes the concentrated likelihood continuous in alpha (with kinks where a unit's maximizing tooth switches) but not in the slopes: where a unit's supremum sits on its feasibility floor (an observation whose count equals its ceiling) the intercept cannot retreat, so when another observation's breakpoint crosses that floor the concentrated likelihood drops by about (1-\alpha)^k, a cliff in the slopes on whose edge the maximizer can rest with a non-zero one-sided gradient, which the perturbed restarts of the outer multistart are there to escape; alternative local optima differ by the order of the jumps (a few hundredths of a log-likelihood unit on a 12-unit panel), far inside the bootstrap variability of the estimates.

Short panels: the dispersion parameter alpha and the fixed effects are subject to the incidental-parameters bias of nonlinear fixed-effects estimation. The bias in alpha is downward and of order 1/T, where T is the per-unit number of observations. It is small for T at least about 30 and should be treated cautiously for short panels. The covariate coefficients carry a smaller bias of the same order, which can matter in very short panels (a few positive counts per unit, as in a zero-truncated intensity). bias_correct = "jackknife" removes the leading 1/T term by the split-panel (half-panel) jackknife: the model is refit on the first and second temporal halves of every unit's series (each unit's rows ordered by time, which the correction requires) and the corrected estimate is 2 * full - mean(halves). In the package's fixed-effects bias Monte Carlo (150 units, alpha = 0.5) the alpha bias at T = 6/10/20/40 falls from -0.122/-0.071/-0.036/-0.018 (T times the bias is -0.71 to -0.73 throughout) to +0.008/+0.004/+0.000/-0.001. On the floor support the bias is larger (-0.153/-0.100/-0.057/-0.033), and the correction leaves 23 to 30 percent of it. With correction on, coefficients and alpha are the corrected estimates (the maximum-likelihood values are kept in ⁠$uncorrected⁠), the unit effects and fitted values are re-concentrated at the corrected parameters, and logLik/AIC continue to refer to the maximum-likelihood fit. The covariate coefficients are corrected too: unlike the Poisson's, whose slopes carry no first-order bias, the CPB's can carry one (see Short panels above). The correction leaves the estimator's asymptotic variance unchanged (Dhaene and Jochmans 2015) but adds variance in short panels. In the package's simulations (40 units, alpha = 0.5, a slope with no bias to remove) the corrected slope's standard deviation was 19 percent above the uncorrected slope's at T = 6 and 13 percent at T = 12, while alpha's bias fell from -0.13 and -0.06 to +0.011 and +0.004; where the slopes' own bias is small next to their standard error, the uncorrected fit (bias_correct = "none") estimates them more precisely. With se = "bootstrap" each replicate is corrected the same way (its two half-panel refits; the validity gate below is a check on the data, made once), so the standard errors include that variance and the bootstrap takes about three times as long. In the same simulations the corrected slope's 95% interval covered 0.91 at T = 6 and 0.94 at T = 12 (0.85 and 0.89 with uncorrected replicates). On the floor support the replicates are not corrected, and the interval is centered on the corrected estimate.

Validity gate: the split-panel identity requires the two half-panels to estimate the same pseudo-true parameter (time-homogeneity; Dhaene & Jochmans 2015, Review of Economic Studies). On trending or time-heterogeneous panels the correction is invalid, so the function refuses it – returning the uncorrected maximum-likelihood fit with a warning naming the failed check – whenever the full fit's residuals trend with time (correlation with the within-unit time index, p < 0.01), the two half-panel dispersion estimates disagree by more than 2.5 times the noise of five placebo splits of the units (and by more than 0.05 of the estimate's scale), or the corrected dispersion leaves its feasible space. A refusal is diagnostic information about the panel, not an error. The gate is deliberately powered over sized: in the package's calibration Monte Carlo it refuses 9 to 13 percent of genuinely time-homogeneous panels (a conservative nuisance; the returned fit is exactly the ordinary maximum-likelihood estimate) while catching 98.5 percent of dispersion regime changes and 100 percent of smooth unmodeled trends – the cases where an uncaught correction would be silently wrong.

One set of fixed effects is concentrated out. For two-way (e.g. unit and time) fixed effects, add factor(time) to formula (entered as dummies), or use count_reg(), which supports two-way fixed effects natively.

Value

An object of class "cpb_fe" with coefficients, alpha, loglik, fe (the estimated unit intercepts), fitted.values, and the implied per-observation ceiling.

See Also

cpb(), ud_screen()

Examples

set.seed(1)
d <- do.call(rbind, lapply(1:25, function(i) {
  x <- rnorm(12); lam <- exp(rnorm(1, 0, 0.5) + 0.5 * x)
  N <- pmax(round(lam / 0.5), 1)
  data.frame(unit = i, x = x, y = rbinom(12, N, 0.5))
}))
fit <- cpb_fe(y ~ x, data = d, fe = "unit")
fit

Cross-validated proper scores for a count model

Description

K-fold cross-validated log score and RPS: the data are split into k folds, the model is refit on each training set via fitfun, and its predicted distribution is scored on the held-out fold. Because the score is out of sample it penalizes over-parameterization, so it is the appropriate criterion for comparing models that differ in the number of parameters (for example a hurdle with per-unit participation fixed effects against a zero-inflated model). Per-observation held-out scores are returned so two models fit on the same folds can be compared with a paired, cluster-robust test. Each fold is scored as score() scores newdata: a unit, or a level of any factor, that the training folds lack enters at the reference unit, and a unit a two-part fit holds at a limit is scored at that limit.

Usage

cv_score(
  fitfun,
  data,
  k = 5,
  kmax = NULL,
  folds = NULL,
  cores = 1L,
  weights = NULL,
  offset = NULL
)

Arguments

Value

An object of class "cv_score": a list with the mean held-out logscore and rps, the per-observation vectors logscore_i/rps_i, and the folds used.

See Also

score()

Examples

set.seed(1)
d <- data.frame(y = rcpb(500, lambda = 3, alpha = 0.5))
cv <- cv_score(function(tr) cpb(y ~ 1, tr, truncated = FALSE, se = "none"), d, k = 5)
cv$logscore

Dispersion profile: variance-to-mean ratio against the fitted mean

Description

Bins the observations by the fitted mean of the first model, computes the empirical conditional variance-to-mean ratio in each bin (the mean of the Pearson contributions (y_i - \hat\mu_i)^2/\hat\mu_i), and sets it against the ratio each fitted family implies at those means. The families differ in how their dispersion moves with the mean: the CPB and the Katz family hold it essentially constant (at alpha and delta; the implied curves are the exact moments of the fitted pmf, so the renormalized CPB's ratio rises above alpha where the implied ceiling is small: by up to about a third at alpha = 0.3 where \lambda/(1-\alpha) lies between one and three, and toward one, the Poisson's value, as the rate falls toward zero; on the floor support it dips below alpha instead), the COM-Poisson, gamma-count, and double Poisson approach a constant only as the mean grows and rise toward the Poisson at small means, and the negative binomial's ratio rises linearly; the hyper-Poisson's cannot fall below 1 - 1/\mu (1 - \mu below a mean of one), so it approaches one as the mean grows, and the EWP's follows its two tail weights. The profile therefore shows which variance function the data follow (a fitted distribution does not by itself identify the process behind it), where a family's fails, and whether the observed underdispersion is confined to a range of means. Zero-truncated fits are profiled on their conditional moments. With frequency weights (those of the first model) the bins hold equal weight rather than equal numbers of rows, each row kept whole, every bin statistic is a weighted mean, and n is the bin's weight total, so a weighted fit is profiled as the data it stands for.

Usage

dispersion_profile(
  object,
  ...,
  bins = 10,
  plot = TRUE,
  ylim = NULL,
  main = "Dispersion profile",
  xlab = "Fitted mean",
  ylab = "Conditional Var / Mean"
)

## S3 method for class 'dispersion_profile'
plot(
  x,
  ylim = NULL,
  main = "Dispersion profile",
  xlab = "Fitted mean",
  ylab = "Conditional Var / Mean",
  ...
)

Arguments

Value

A data frame of class "dispersion_profile" with one row per bin: bin, n, mean_fitted, mean_y, ratio_empirical, and one ⁠ratio_<model>⁠ column per fitted model.

See Also

dispersion_test(), compare_models(), rootogram()

Examples

set.seed(4); n <- 400; x <- rnorm(n)
d <- data.frame(y = rgammacount(n, exp(0.8 + 0.6 * x), alpha = 2.5), x = x)
m_gc <- count_reg(y ~ x, d, family = "gammacount")
m_nb <- count_reg(y ~ x, d, family = "negbin")
m_cp <- cpb(y ~ x, d, truncated = FALSE, se = "none")
dispersion_profile(gammacount = m_gc, negbin = m_nb, cpb = m_cp)

Regression-adjusted test of equidispersion

Description

Tests whether a fitted count model's dispersion departs from the Poisson, conditional on the covariates. For a fit from a family with a dispersion parameter (method = "lr", the default), the test is the likelihood-ratio test of that parameter against its Poisson value, with the Poisson refit on the same design (fixed effects included); its p-value comes from the asymptotic distribution or from a parametric bootstrap under the fitted Poisson, by the rule in Details. For a Poisson fit, method = "auxiliary" runs the Cameron–Trivedi (1990) auxiliary regression, which needs no alternative family.

Usage

dispersion_test(
  object,
  alternative = c("two.sided", "under", "over"),
  method = c("lr", "auxiliary"),
  B = NULL,
  cores = 1L
)

Arguments

Details

The Poisson is nested in every family at one value of the dispersion parameter: alpha = 1 (CPB), delta = 1 (GEC), nu = 1 (both COM-Poisson parameterizations), lambda = 0 (generalized Poisson), alpha = 1 (gamma-count), theta = 1 (double Poisson), gamma = 1 (hyper-Poisson), beta1 = beta2 = 0 or beta = 0 (the EWP), and 1/theta = 0 (negative binomial). When that value is interior to the parameter space the likelihood-ratio statistic is asymptotically \chi^2_1 under the null; a directional alternative ("under" or "over") uses the signed root r = \pm\sqrt{LR}, positive when the estimate lies on the underdispersed side of the Poisson value, which is asymptotically standard normal. When the Poisson value is on the boundary (the CPB, whose alpha lives in (0, 1); the negative binomial, whose overdispersion parameter is non-negative) the asymptotic null distribution is the \tfrac12\chi^2_0 + \tfrac12\chi^2_1 mixture (Self and Liang 1987; Andrews 2001) and the test is one-sided by construction ("under" for the CPB, "over" for the negative binomial). The CPB qualifies although its support moves with its parameters: near alpha = 1 its ceiling lambda / (1 - alpha) diverges, the log-likelihood is regular on that side, and its score there is the classical dispersion score -[(y - \lambda)^2 - y] / (2\lambda); the zero-truncated model reaches the same mixture through its efficient score. Because the Poisson is the CPB's limit as alpha approaches 1, the statistic cannot be negative in principle; when the fit reaches that limit (King's CPB, whose likelihood approaches the Poisson's continuously) or the floor support's max.support guard stops alpha short of it, and the fitted log-likelihood lands just below the Poisson's, the statistic takes its boundary value 0. Any other negative statistic means the optimizer fell short, and no p-value is reported.

Calibration. Where the Poisson value is on the boundary of the family's parameter space (the CPB, the negative binomial), or the fit has unit fixed effects, the p-value is by default a parametric bootstrap with 199 replicates, which does not depend on the direction in which the asymptotic reference errs. For the CPB it errs toward rejection: in the package's simulations from the Poisson the asymptotic boundary mixture rejected 7 to 10 percent of 5% tests with 60 to 200 observations, and 54 to 92 percent with unit fixed effects, whose intercepts bias the dispersion estimate toward underdispersion by an amount of order 1/T. The same rule bootstraps the negative binomial's boundary test. Where the Poisson value is interior, the asymptotic distribution is used unless the estimated mean parameters shift it too far: to first order they move the signed root toward underdispersion by p/\sqrt{2n} standard deviations for p mean parameters on n observations (the weight total), the bias Dean and Lawless (1989) remove from the score test, and the bootstrap takes over when that shift pushes the first-order size of a 5% test above 6%. The shift is derived for one tested parameter, so the EWP3's joint test of its two tails takes the bootstrap by default: in the package's simulations from the Poisson its asymptotic \chi^2_2 rejected 6.4 to 7.5 percent of 5% tests with one or two mean parameters on 200 to 500 observations, and 8.0 percent with ten on 400. The bootstrap simulates responses from the Poisson (zero-truncated Poisson) fitted to the same design, offset, weights, and unit effects, refits both models to each, and reports the share of simulated statistics at least as extreme as the observed one, counting the observed one, among the replicates whose two fits both succeeded. Simulating at the Poisson value itself keeps the bootstrap valid at the boundary. B = 0 requests the asymptotic distribution (with unit fixed effects it returns the statistic without a p-value), and a positive B sets the number of bootstrap replicates. Every replicate refits both models, so the bootstrap takes about B times as long as the original fit; cores spreads the replicates over worker processes.

Value

An object of class c("dispersion_test", "htest") with the statistic, the p-value, the estimated dispersion parameter and its Poisson value, the two log-likelihoods, and the alternative; for method = "lr" also calibration ("asymptotic" or "parametric bootstrap"), B and B_ok (replicates requested and completed), boot (the simulated statistics), first_order_size (the asymptotic 5% test's first-order size), asymptotic_ok (whether the calibration rule admits the asymptotic distribution), and note where a statistic or p-value needs one.

References

Andrews, D. W. K. (2001). Testing when a parameter is on the boundary of the maintained hypothesis. Econometrica, 69(3), 683-734. Cameron, A. C. and Trivedi, P. K. (1990). Regression-based tests for overdispersion in the Poisson model. Journal of Econometrics, 46(3), 347-364. Dean, C. and Lawless, J. F. (1989). Tests for detecting overdispersion in Poisson regression models. Journal of the American Statistical Association, 84(406), 467-472. Self, S. G. and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82(398), 605-610.

See Also

count_reg(), dispersion_profile(), zi_test()

Examples

set.seed(2); n <- 400; x <- rnorm(n)
d <- data.frame(y = rgammacount(n, exp(1 + 0.4 * x), alpha = 2), x = x)
dispersion_test(count_reg(y ~ x, d, family = "gammacount"), alternative = "under")
dispersion_test(count_reg(y ~ x, d, family = "poisson"), method = "auxiliary")

Double Poisson distribution functions

Description

Density, distribution, quantile, and random generation for Efron's (1986) double Poisson distribution with ___location mu and dispersion theta, normalized by the exact sum rather than Efron's closed-form approximation. theta > 1 gives underdispersion, theta = 1 is the Poisson, theta < 1 overdispersion; the variance-to-mean ratio is approximately 1/theta, and the mean is close to (but not exactly) mu. The exact mean and variance are what count_reg(family = "doublepois") reports through predict() and fitted().

Usage

ddoublepois(x, mu, theta, log = FALSE)

pdoublepois(q, mu, theta, lower.tail = TRUE, log.p = FALSE)

qdoublepois(p, mu, theta, lower.tail = TRUE, log.p = FALSE)

rdoublepois(n, mu, theta)

Arguments

Value

ddoublepois a density, pdoublepois a CDF, qdoublepois a quantile, rdoublepois a numeric vector of count draws.

References

Efron, B. (1986). Double exponential families and their use in generalized linear regression. Journal of the American Statistical Association, 81(395), 709-721.

See Also

count_reg()

Examples

ddoublepois(0:5, mu = 3, theta = 2)
sum(ddoublepois(0:60, mu = 3, theta = 2))

Exponentially weighted Poisson distribution functions

Description

Density, distribution, quantile, and random generation for the exponentially weighted Poisson of Ridout and Besbeas (2004), in its three-parameter form (EWP3) and, with beta2 = beta1, its two-parameter form (EWP2): a Poisson with centrality lambda reweighted by a two-sided exponential in the distance from it,

P(Y = y) \propto w_y e^{-\lambda} \lambda^y / y!, \quad w_y = e^{-\beta_1 (\lambda - y)} \ (y \le \lambda), \quad w_y = e^{-\beta_2 (y - \lambda)} \ (y > \lambda).

beta1 and beta2 greater than zero thin the lower and upper tails (underdispersion), values below zero thicken them, and beta1 = beta2 = 0 is the Poisson. lambda is not the mean; the normalizer and the exact mean and variance are sums over the support (no truncation of the normalizing sum at a fixed count), and the exact mean is what count_reg(family = "ewp3") and count_reg(family = "ewp2") report through predict() and fitted().

Usage

dewp(x, lambda, beta1, beta2 = beta1, log = FALSE)

pewp(q, lambda, beta1, beta2 = beta1, lower.tail = TRUE, log.p = FALSE)

qewp(p, lambda, beta1, beta2 = beta1, lower.tail = TRUE, log.p = FALSE)

rewp(n, lambda, beta1, beta2 = beta1)

Arguments

Value

dewp a density, pewp a CDF, qewp a quantile, rewp a numeric vector of count draws.

References

Ridout, M. S. and Besbeas, P. (2004). An empirical model for underdispersed count data. Statistical Modelling, 4(1), 77-89.

See Also

count_reg()

Examples

dewp(0:6, lambda = 3, beta1 = 1, beta2 = 1)
var(rewp(2000, lambda = 3, beta1 = 1, beta2 = 1))

First difference for a CPB fit

Description

The effect on a quantity of interest of moving one covariate from one value to another, holding the other covariates at their means (or at the values at gives them).

Usage

first_difference(object, ...)

## S3 method for class 'cpb'
first_difference(
  object,
  variable,
  from,
  to,
  quantity = c("mean", "ceiling", "prob"),
  y = NULL,
  level = 0.95,
  at = NULL,
  ...
)

Arguments

Details

Every first_difference() method in the package takes its interval by one rule. When the fit carries bootstrap draws (se = "bootstrap" at fit time), the quantity is computed at every draw, jointly over the coefficients and the dispersion parameter (and, for a fixed-effects fit, the reference unit's effect), and the percentile interval is taken; the draws of a fixed-effects fit are first recentered on the reported estimate, because its replicates share the incidental-parameters bias. Otherwise the interval is the delta method's over the fit's full covariance, the dispersion parameter included (the matched count families' analytic, robust or clustered covariance; a hurdle's two margins, which are independent). A fit with neither gets the point estimate with method = "none".

Value

A "ud_fd" data frame – the package-wide first-difference contract (columns component, from, to, diff, lower, upper, method) – with one row for the requested quantity and a bootstrap percentile interval on the difference (method = "bootstrap (stored)"). Every first_difference() method in the package returns this same shape.

Examples

set.seed(1); x <- rnorm(250)
N <- pmax(round(exp(1.6 + 0.5 * x) / 0.5), 1); y <- rbinom(250, N, 0.5)
fit <- cpb(y ~ x, data = data.frame(y = y, x = x)[y > 0, ], se = "bootstrap", B = 100)
first_difference(fit, "x", from = -1, to = 1, quantity = "mean")

First differences for the matched count families

Description

The discrete-change effect on the expected count E(Y) as variable moves from from to to, holding the other covariates at their sample means (or at the values at gives them). The interval follows the package's one rule (see first_difference()): percentile intervals from the fit's stored bootstrap draws for the GEC, fixed-effects CPB and zero-inflated GEC fits, the delta method over the full covariance, the dispersion parameter included, for the matched count families and the hurdles' margins. For count_reg this is a single number; for the two-part models it is reported on each margin: the binary stage's probability (participation, or the structural-zero probability), the count stage's mean, and the marginal mean, their product, each at from and to. The first two are different quantities and their changes do not add up to the marginal change, which is \Delta(pm) = \Delta p\, m_0 + p_0 \Delta m + \Delta p\, \Delta m.

Usage

## S3 method for class 'gec'
first_difference(object, variable, from, to, level = 0.95, at = NULL, ...)

## S3 method for class 'gec_fe'
first_difference(object, variable, from, to, level = 0.95, at = NULL, ...)

## S3 method for class 'cpb_fe'
first_difference(object, variable, from, to, level = 0.95, at = NULL, ...)

## S3 method for class 'hurdle_gec'
first_difference(
  object,
  variable,
  from,
  to,
  level = 0.95,
  stage = c("both", "participation", "intensity", "zero", "count"),
  at = NULL,
  ...
)

## S3 method for class 'zi_gec'
first_difference(
  object,
  variable,
  from,
  to,
  level = 0.95,
  stage = c("both", "participation", "intensity", "zero", "count"),
  at = NULL,
  ...
)

## S3 method for class 'count_reg'
first_difference(object, variable, from, to, level = 0.95, at = NULL, ...)

## S3 method for class 'hurdle_count'
first_difference(
  object,
  variable,
  from,
  to,
  level = 0.95,
  stage = c("both", "participation", "intensity", "zero", "count"),
  at = NULL,
  ...
)

## S3 method for class 'zi_count'
first_difference(
  object,
  variable,
  from,
  to,
  level = 0.95,
  stage = c("both", "participation", "intensity", "zero", "count"),
  at = NULL,
  ...
)

Arguments

Value

A "ud_fd" data frame – the package-wide first-difference contract shared by every first_difference() method: columns component, from, to, diff, lower, upper, method. Single-equation fits return one row (component = "mean"); two-part fits return one row per margin level (binary-stage probability, count-stage mean, marginal mean). method records the uncertainty source per row ("delta", a bootstrap label, or "none" with NA bounds).

Extensive/intensive first-difference decomposition for a hurdle- or ZI-CPB

Description

The effect of moving variable from one value to another on the participation (extensive) margin, the intensity (intensive) margin, and the marginal expected count, holding other covariates at their means (or at the values at gives them). The two margins are reported side by side; their changes do not add up to the marginal change (see first_difference.count). With B = 0 the intervals follow the package's one rule (see first_difference()): the stored bootstrap draws of a zi_cpb fit, and for a hurdle the delta method over its two margins' covariances (the participation model's, and the intensity's stored bootstrap draws). B > 0 refits the whole model on resampled data instead.

Usage

## S3 method for class 'hurdle_cpb'
first_difference(
  object,
  variable,
  from,
  to,
  B = 0,
  level = 0.95,
  data = NULL,
  stage = c("both", "participation", "intensity", "zero", "count"),
  cores = 1L,
  at = NULL,
  ...
)

## S3 method for class 'zi_cpb'
first_difference(
  object,
  variable,
  from,
  to,
  B = 0,
  level = 0.95,
  data = NULL,
  stage = c("both", "participation", "intensity", "zero", "count"),
  cores = 1L,
  at = NULL,
  ...
)

Arguments

Value

A "ud_fd" data frame (the package-wide first-difference contract): columns component, from, to, diff, lower, upper, method, one row per margin level.

Fitted means for the two-part and fixed-effects classes

Description

Completes stats::fitted() across the family, returning the same quantity as predict(object, type = "response") on the estimation data: the mean for the single-equation fits and the marginal expected count (zeros included) for the two-part fits. gec_fe inherits fitted.gec.

Usage

## S3 method for class 'cpb_fe'
fitted(object, ...)

## S3 method for class 'hurdle_cpb'
fitted(object, ...)

## S3 method for class 'hurdle_gec'
fitted(object, ...)

## S3 method for class 'hurdle_count'
fitted(object, ...)

## S3 method for class 'zi_cpb'
fitted(object, ...)

## S3 method for class 'zi_gec'
fitted(object, ...)

## S3 method for class 'zi_count'
fitted(object, ...)

Arguments

Value

A numeric vector, one element per estimation observation.

Gamma-count distribution functions

Description

Density, distribution, quantile, and random generation for Winkelmann's (1995) gamma-count distribution: the number of events in a unit interval when the waiting times between events are independent Gamma variables with shape alpha and rate alpha * mu, so that mu is the long-run event rate. alpha > 1 (waiting times more regular than exponential) gives underdispersion, alpha = 1 is the Poisson, alpha < 1 overdispersion; the variance-to-mean ratio is approximately 1/alpha. The exact mean and variance are finite sums of incomplete-gamma terms and are what count_reg(family = "gammacount") reports through predict() and fitted().

Usage

dgammacount(x, mu, alpha, log = FALSE)

pgammacount(q, mu, alpha, lower.tail = TRUE, log.p = FALSE)

qgammacount(p, mu, alpha, lower.tail = TRUE, log.p = FALSE)

rgammacount(n, mu, alpha)

Arguments

Value

dgammacount a density, pgammacount a CDF, qgammacount a quantile, rgammacount a numeric vector of count draws.

References

Winkelmann, R. (1995). Duration dependence and dispersion in count-data models. Journal of Business & Economic Statistics, 13(4), 467-474.

See Also

count_reg()

Examples

dgammacount(0:5, mu = 3, alpha = 2)
var(rgammacount(2000, mu = 3, alpha = 2)) / 3

Generalized event count (Katz-family) regression

Description

Fits King's generalized event count model, a Katz-family count regression whose single dispersion parameter delta (the variance-to-mean ratio) is estimated freely and spans underdispersion (delta < 1, a finite-support member, the continuous parameter binomial of cpb() at alpha = delta), equidispersion (delta = 1, Poisson), and overdispersion (delta > 1, the negative binomial). Unlike cpb(), which fixes the direction of dispersion to under, gec() lets the data choose; wherever its estimate of delta is below one the two return the same fit. The Katz recursion delivers exact first and second moments, summed over the support (finite for delta < 1, where they are the correction of Winkelmann, Signorino and King), and the likelihood is evaluated in C++. For delta < 1 the support is finite and the renormalized distribution's mean and variance equal ⁠exp(x'b)⁠ and ⁠delta * exp(x'b)⁠ only when ⁠exp(x'b)/(1 - delta)⁠ is an integer, so ⁠exp(x'b)⁠ is the rate parameter of the recursion, the fitted mean is computed exactly from the pmf, and delta is the Katz dispersion parameter (the variance-to-mean ratio on an unbounded support). The optimizer works on covariates scaled to unit standard deviation and maps the coefficients back, so the fit does not depend on the covariates' units.

Usage

gec(
  formula,
  data,
  truncated = FALSE,
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  offset = NULL,
  weights = NULL,
  cores = 1L,
  max.support = NULL,
  maxit = 20000,
  reltol = 1e-08
)

Arguments

Value

An object of class "gec" with coefficients, delta (the estimated dispersion), loglik, bootstrap standard errors, and bookkeeping.

See Also

cpb(), compare_dispersion(), dispersion_test()

Examples

set.seed(1); x <- rnorm(400)
y <- rpois(400, exp(1 + 0.5 * x))
gec(y ~ x, data = data.frame(y = y, x = x), se = "none")

Generalized event count (Katz) distribution functions

Description

Density, distribution, quantile, and random generation for the King (1989) generalized event count model with rate lambda and Katz dispersion delta (⁠< 1⁠ underdispersed, ⁠= 1⁠ Poisson, ⁠> 1⁠ overdispersed). On an unbounded support (delta >= 1) the mean is lambda and the variance-to-mean ratio is delta exactly; for delta < 1 the support is finite and the renormalized distribution's mean and variance equal those values only when lambda/(1 - delta) is an integer (the exact moments are finite sums of the pmf, as count_reg()'s siblings report them). Consistent with the estimator gec().

Usage

dgec(x, lambda, delta, max.support = 500, log = FALSE)

pgec(q, lambda, delta, max.support = 500, lower.tail = TRUE, log.p = FALSE)

qgec(p, lambda, delta, max.support = 500, lower.tail = TRUE, log.p = FALSE)

rgec(n, lambda, delta, max.support = 500)

Arguments

Value

dgec a density, pgec a CDF, qgec a quantile, rgec a numeric vector of count draws.

See Also

gec()

Examples

dgec(0:5, lambda = 3, delta = 0.7)
var(rgec(2000, lambda = 3, delta = 0.7)) / mean(rgec(2000, lambda = 3, delta = 0.7))

GEC (Katz-family) regression with high-dimensional unit fixed effects

Description

Fits gec() with a full set of unit fixed effects, concentrating (profiling) out the unit intercepts by a one-dimensional inner maximization per unit, so the outer optimizer handles only the covariate coefficients and the free dispersion parameter delta. This is the cpb_fe() concentration generalized to the whole Katz family: the panel can be under-, equi-, or overdispersed and the direction is estimated, not presumed. For delta < 1 the Katz support is finite, but the likelihood is continuous across an integer ceiling (the entering support point's mass grows from zero) and only kinks there, so each unit's objective is unimodal in its intercept and is maximized by golden section on a bracket around the unit's Poisson intercept, expanded while the maximum sits at an edge; a maximum at a kink is tracked in the analytic gradient of the concentrated likelihood, which the outer BFGS uses. The feasibility floor is exact: every count must lie in ⁠0..ceiling(mu/(1-delta))⁠. Runtime: about one and a half times that of cpb_fe() on the same panel (the Katz recursion runs the whole support): on the positive counts of the bundled peacekeeping panel, 2,602 rows in 146 units, gec_fe() takes about two and a half minutes and cpb_fe() about one and a half, and the two return the same fit.

Usage

gec_fe(
  formula,
  data,
  fe,
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  offset = NULL,
  weights = NULL,
  cores = 1L,
  max.support = NULL,
  inner_it = 30L,
  maxit = 3000L,
  reltol = 1e-07,
  bias_correct = c("none", "jackknife"),
  time = NULL
)

Arguments

Details

Limitation: unlike cpb_fe(), gec_fe() has no truncated argument – a concentrated zero-truncated GEC is not currently implemented. For a zero-truncated GEC with fixed effects, enter the unit factor as dummies in gec()'s formula (feasible for moderate unit counts).

Short panels: delta and the fixed effects carry the incidental-parameters bias of nonlinear fixed-effects estimation, of order 1/T; treat delta cautiously below about T = 30. The covariate coefficients carry a smaller bias of the same order. bias_correct = "jackknife" removes the leading 1/T term by the split-panel jackknife exactly as in cpb_fe(): refit on each unit's temporal halves and report 2 * full - mean(halves), with the unit effects and fitted values re-concentrated at the corrected parameters and logLik/AIC kept at the maximum-likelihood fit (uncorrected estimates in ⁠$uncorrected⁠). The same time-homogeneity validity gate applies: on panels where the two halves do not estimate a common parameter the correction is refused with a warning and the maximum-likelihood fit is returned (see cpb_fe(), Details). As there, the covariate coefficients are corrected too, and with se = "bootstrap" each replicate is corrected the same way, so the standard errors include the correction's own variance.

Value

An object of class c("gec_fe", "gec").

See Also

gec(), cpb_fe()

Examples

set.seed(5)
d <- do.call(rbind, lapply(1:30, function(i) {
  x <- rnorm(10); N <- pmax(round(exp(1 + rnorm(1, 0, 0.4) + 0.3 * x) / 0.5), 1)
  data.frame(unit = i, x = x, y = rbinom(10, N, 0.5))
}))
gec_fe(y ~ x, data = d, fe = "unit")   # delta ~ 0.5: underdispersed

Generalized Poisson distribution functions

Description

Density, distribution, quantile, and random generation for the Consul-Jain generalized Poisson distribution in its constant-dispersion form: mean mu and dispersion lambda in (-1, 1), with P(Y = y) = \theta(\theta + \lambda y)^{y-1} e^{-\theta - \lambda y}/y! and \theta = \mu(1 - \lambda), so that the variance-to-mean ratio is 1/(1-\lambda)^2. lambda < 0 gives underdispersion, lambda = 0 the Poisson, lambda > 0 overdispersion. For lambda < 0 the support is finite (the terms are positive only while \theta + \lambda y > 0) and the pmf is renormalized on that support, so it is a proper distribution and the exact mean and variance (which then differ from mu and mu/(1-lambda)^2 only by the negligible mass beyond the support) are what count_reg(family = "genpois") reports through predict() and fitted().

Usage

dgenpois(x, mu, lambda, log = FALSE)

pgenpois(q, mu, lambda, lower.tail = TRUE, log.p = FALSE)

qgenpois(p, mu, lambda, lower.tail = TRUE, log.p = FALSE)

rgenpois(n, mu, lambda)

Arguments

Value

dgenpois a density, pgenpois a CDF, qgenpois a quantile, rgenpois a numeric vector of count draws.

References

Consul, P. C. and Jain, G. C. (1973). A generalization of the Poisson distribution. Technometrics, 15(4), 791-799. Consul, P. C. and Famoye, F. (2006). Lagrangian Probability Distributions. Birkhauser.

See Also

count_reg()

Examples

dgenpois(0:5, mu = 3, lambda = -0.3)
var(rgenpois(2000, mu = 3, lambda = -0.3)) / 3   # about 1/(1.3)^2

Glance at a fitted model (broom method)

Description

One row of fit statistics for every model class of the package: logLik, AIC, BIC, df and nobs (the weight total of a weighted fit), the dispersion parameter as dispersion with its name (dispersion_parameter) and the direction its estimate points (dispersion_direction; the families' parameters point different ways, so the CPB's alpha below 1 and the gamma-count's alpha above 1 are both underdispersion), coef_scale, the scale of the count equation's coefficients (see tidy.cpb()), and the class's own columns (alpha for the CPB family, family for the matched count families, n.units for the fixed-effects fits).

Usage

## S3 method for class 'cpb'
glance(x, ...)

## S3 method for class 'cpb_fe'
glance(x, ...)

## S3 method for class 'hurdle_cpb'
glance(x, ...)

## S3 method for class 'zi_cpb'
glance(x, ...)

## S3 method for class 'gec'
glance(x, ...)

## S3 method for class 'hurdle_gec'
glance(x, ...)

## S3 method for class 'zi_gec'
glance(x, ...)

## S3 method for class 'count_reg'
glance(x, ...)

## S3 method for class 'zi_count'
glance(x, ...)

## S3 method for class 'hurdle_count'
glance(x, ...)

Arguments

Value

A one-row data frame of fit statistics.

Hurdle count regression for the matched count families

Description

Fits a participation model joined to a zero-truncated intensity from any of the package's count families (see count_reg(), Families), the count analogue of hurdle_cpb(). Returned as a "hurdle_count" object that compare_models() and score() accept, so a hurdle-CPB and a hurdle gamma-count can be compared on one footing.

Usage

hurdle_count(
  formula,
  data,
  family = c("poisson", "negbin", "compois", "mpcmp", "genpois", "gammacount",
    "doublepois", "hyperpois", "ewp3", "ewp2"),
  participation = NULL,
  fe = NULL,
  part_fe = NULL,
  link = c("logit", "probit", "cloglog"),
  offset = NULL,
  weights = NULL,
  se = c("analytic", "robust", "cluster", "none"),
  cluster = NULL
)

Arguments

Value

An object of class "hurdle_count".

See Also

hurdle_cpb(), zi_count(), compare_models()

Examples

set.seed(1); n <- 300; x <- rnorm(n); z <- rnorm(n)
y <- ifelse(rbinom(n, 1, plogis(0.4 + 0.8 * z)) == 1, rpois(n, exp(1 + 0.3 * x)) + 1L, 0L)
hurdle_count(y ~ x, data.frame(y = y, x = x, z = z), family = "poisson",
             participation = ~ z, se = "none")

Hurdle continuous parameter binomial regression

Description

Fits a participation (hurdle) model joined to a zero-truncated CPB for the positive counts. Because the hurdle log-likelihood factorizes, the two parts are fit separately: a logistic regression of participation on all units, and a zero-truncated CPB (cpb()) on the positive counts. The result is the natural model for a bounded, underdispersed participation process — most units at zero, participants carrying a tight count.

Usage

hurdle_cpb(
  formula,
  data,
  participation = NULL,
  fe = NULL,
  part_fe = NULL,
  link = c("logit", "probit", "cloglog"),
  offset = NULL,
  cluster = NULL,
  se = c("none", "bootstrap"),
  B = 500,
  weights = NULL,
  cores = 1L,
  ...
)

Arguments

Value

An object of class "hurdle_cpb" with elements participation (a glm), intensity (a cpb or cpb_fe), and bookkeeping.

Examples

set.seed(1)
n <- 800; x <- rnorm(n); z <- rnorm(n)
y <- rhurdle_cpb(n, lambda = exp(1.3 + 0.5 * x), alpha = 0.5,
                 p = plogis(-0.2 + 0.8 * z))
fit <- hurdle_cpb(y ~ x, data = data.frame(y = y, x = x, z = z),
                  participation = ~ z)
fit

Hurdle GEC (Katz-family) regression

Description

Joins a participation (hurdle) logit to a zero-truncated GEC intensity on the positive counts. Unlike hurdle_cpb(), whose intensity is fixed to the underdispersed CPB, the GEC intensity estimates its own dispersion direction.

Usage

hurdle_gec(
  formula,
  data,
  participation = NULL,
  part_fe = NULL,
  link = c("logit", "probit", "cloglog"),
  offset = NULL,
  cluster = NULL,
  se = c("none", "bootstrap"),
  B = 500,
  weights = NULL,
  cores = 1L,
  max.support = NULL
)

Arguments

Value

An object of class "hurdle_gec".

Limitation

Unlike hurdle_cpb(), hurdle_gec() has no fe argument for a concentrated fixed-effects intensity (no concentrated zero-truncated GEC is implemented). For within-unit intensities, enter the unit factor as dummies in formula (feasible for moderate unit counts), or use hurdle_cpb().

See Also

hurdle_cpb(), gec(), zi_gec()

Examples

set.seed(3); n <- 300; x <- rnorm(n); z <- rnorm(n)
y <- ifelse(rbinom(n, 1, plogis(0.3 + 0.8 * z)) == 1, rpois(n, exp(1 + 0.3 * x)) + 1L, 0L)
hurdle_gec(y ~ x, data.frame(y = y, x = x, z = z), participation = ~ z)

Hyper-Poisson distribution functions

Description

Density, distribution, quantile, and random generation for the hyper-Poisson distribution of Bardwell and Crow (1964) in the mean parameterization of Saez-Castillo and Conde-Sanchez (2013): mean mu and dispersion gamma > 0, with P(Y = y) = \lambda^y / ((\gamma)_y\, {}_1F_1(1; \gamma; \lambda)), (\gamma)_y = \Gamma(\gamma + y)/\Gamma(\gamma), and the rate \lambda the unique solution of E(Y) = \mu. gamma < 1 gives underdispersion, gamma = 1 the Poisson, gamma > 1 overdispersion. The underdispersion is bounded: as gamma -> 0 the distribution tends to 1 + Poisson(mu - 1) for mu > 1 and to a Bernoulli(mu) below, so the variance is at least mu - 1 (mu >= 1) or mu (1 - mu) (mu < 1), and Var/Mean approaches one as the mean grows (0.5 at a mean of 2, 0.9 at 10). The normalizer and the moments are exact sums over the support.

Usage

dhyperpois(x, mu, gamma, log = FALSE)

phyperpois(q, mu, gamma, lower.tail = TRUE, log.p = FALSE)

qhyperpois(p, mu, gamma, lower.tail = TRUE, log.p = FALSE)

rhyperpois(n, mu, gamma)

Arguments

Value

dhyperpois a density, phyperpois a CDF, qhyperpois a quantile, rhyperpois a numeric vector of count draws.

References

Bardwell, G. E. and Crow, E. L. (1964). A two-parameter family of hyper-Poisson distributions. Journal of the American Statistical Association, 59(305), 133-141. Sáez-Castillo, A. J. and Conde-Sánchez, A. (2013). A hyper-Poisson regression model for overdispersed and underdispersed count data. Computational Statistics & Data Analysis, 61, 148-157.

See Also

count_reg()

Examples

dhyperpois(0:5, mu = 3, gamma = 0.5)
var(rhyperpois(2000, mu = 3, gamma = 0.5)) / 3

Implied ceiling, and the ceiling at the limits of alpha's interval

Description

Returns the observation- (or profile-) specific ceiling lambda/(1-alpha), and the ceiling at the fitted rate with alpha at the limits of its interval (see alpha_confint(); calibrated when the fit went through calibrate_alpha()). Those two columns hold the rate at its estimate, so they show how far the dispersion's uncertainty alone moves the ceiling and are not a confidence interval for it; first_difference() with quantity = "ceiling" propagates the coefficients' uncertainty too. The ceiling is the CPB's continuous number of trials. No count above its integer ceiling has positive probability on King's support (the default), and none above its integer floor on the floor support (support = "floor"); that largest feasible count is returned as max_count. Whether held-out counts respect the bound is tested by ceiling_exceedance().

Usage

implied_ceiling(object, ...)

## S3 method for class 'cpb'
implied_ceiling(object, newdata = NULL, level = 0.95, ...)

## S3 method for class 'cpb_fe'
implied_ceiling(object, newdata = NULL, level = 0.95, ...)

Arguments

Details

The ceiling is the fitted distribution's support end. With the rate and alpha it is fixed by the CPB's variance-mean link, so it is a bound on the counts only when the CPB is the data's process. On counts with a common hard bound from another process (binomial counts out of six trials, say), the fitted ceiling follows the mean and the dispersion instead: in one simulation it ran from 3 where the mean was low to 12 where it was high, below the true bound in the first case and well above it in the second.

Value

A data frame with lambda, ceiling, max_count (the largest count with positive probability), and ceiling_at_alpha_lower / ceiling_at_alpha_upper, the ceiling at the fitted rate with alpha at the lower and upper limits of its interval (for cpb fits, and hurdles with a cpb intensity).

Examples

set.seed(6); x <- rnorm(300)
N <- pmax(round(exp(1.5 + 0.4 * x) / 0.5), 1); y <- rbinom(300, N, 0.5)
fit <- cpb(y ~ x, data.frame(y = y, x = x)[y > 0, ], se = "none")
implied_ceiling(fit, newdata = data.frame(x = c(-1, 0, 1)))

Incidence rate ratios for a CPB fit

Description

The ratios are \exp(\beta), ratios of the CPB's rate parameter \lambda. The CPB's mean equals \lambda only where \lambda/(1-\alpha) is an integer, and a zero-truncated fit (the cpb() default) models E(Y \mid Y \ge 1), so a ratio of expected counts can differ from them; first_difference() reports the exact mean contrast.

Usage

irr(object, ...)

## S3 method for class 'cpb'
irr(object, level = 0.95, ...)

Arguments

Value

A "ud_irr" data frame – the package-wide rate-ratio contract (columns term, equation, ratio, estimate, lower, upper, method) shared by every irr() method; here with bootstrap percentile intervals from the stored draws.

Examples

set.seed(8); x <- rnorm(200)
N <- pmax(round(exp(1.5 + 0.4 * x) / 0.5), 1); y <- rbinom(200, N, 0.5)
fit <- cpb(y ~ x, data.frame(y = y, x = x)[y > 0, ], se = "bootstrap", B = 100)
irr(fit)

Rate and odds ratios for the matched count families

Description

Rate ratios (exp of the count/intensity coefficients, column IRR) and, for the two-part models, odds ratios of the binary stage (column OR), by equation. Element names match the CPB family (intensity for the count component; participation for the hurdle stage, zero for the inflation stage).

Usage

## S3 method for class 'count_reg'
irr(object, level = 0.95, ...)

## S3 method for class 'gec'
irr(object, level = 0.95, ...)

## S3 method for class 'gec_fe'
irr(object, level = 0.95, ...)

## S3 method for class 'cpb_fe'
irr(object, level = 0.95, ...)

## S3 method for class 'hurdle_gec'
irr(object, level = 0.95, ...)

## S3 method for class 'zi_gec'
irr(object, level = 0.95, ...)

## S3 method for class 'hurdle_count'
irr(object, level = 0.95, ...)

## S3 method for class 'zi_count'
irr(object, level = 0.95, ...)

Arguments

Value

A "ud_irr" data frame – the package-wide ratio contract shared by every irr() method: columns term, equation ("count", "participation" for a hurdle's binary stage, or "inflation" for a zero-inflated model's structural-zero stage; the two point in opposite directions, an odds ratio above one raising P(Y > 0) in the first and P(structural zero) in the second), ratio ("IRR", "OR", or the parameter's ratio named below), estimate, lower, upper, and method (the interval source; "none" with NA bounds when no covariance is available). Two-part models stack both equations. The count rows are \exp(\beta), the ratio of the family's natural parameter. That is a ratio of means for the Poisson, the negative binomial, the mean-parameterized COM-Poisson ("mpcmp"), the hyper-Poisson, and the generalized Poisson with \lambda \ge 0 ("IRR"). It is not for the others, whose rows name what it is, with a message: the rate-parameterized COM-Poisson ("compois", "lambda ratio"; its mean moves roughly as the rate to the power 1/\nu, so that \exp(\beta) is roughly the mean ratio to the power \nu), the EWP ("centrality ratio"), the gamma count ("long-run rate ratio"), the double Poisson ("___location ratio") and the generalized Poisson below zero ("mu ratio"). The CPB and the underdispersed GEC keep "IRR", a ratio of their rate, which their mean follows only approximately (the rate renormalized over a finite support). For all of these, first_difference() and ame() report the exact mean contrast.

Mundlak (correlated random effects) device

Description

Augments a data frame with the unit-level means of the time-varying numeric covariates in formula, and returns the augmented formula and data. Fitting any of the package's estimators on the result implements the Mundlak / correlated-random-effects specification: the coefficients on the original covariates recover the within-unit (fixed-effects-consistent) effects, while the coefficients on the unit means capture (and test) the correlation between the covariates and the unit effect. It is a lighter, between-variation-preserving alternative to full fixed effects, and it is model-agnostic — the same device feeds cpb(), gec(), or count_reg().

Usage

mundlak(formula, data, unit, suffix = "_mean")

Arguments

Value

A list with formula (the augmented formula), data (the augmented data frame), and added (the names of the unit-mean columns).

See Also

count_reg(), cpb()

Examples

set.seed(1)
d <- data.frame(y = rpois(200, 3), x = rnorm(200), unit = factor(rep(1:20, each = 10)))
m <- mundlak(y ~ x, d, unit = "unit")
fit <- count_reg(m$formula, m$data, family = "poisson")

UN peacekeeping contributions, state-years

Description

A state-year panel of contributions to United Nations peacekeeping operations, the archetype of a bounded, underdispersed participation count: most state-years contribute to no operation, and the states that contribute carry a tight count held in a narrow band by a turning-over roster. About 42% of state-years are zero, so it is also a natural example for the hurdle-CPB.

Usage

data(peacekeeping)

Format

A data frame with 4,448 rows and 9 variables:

iso3

ISO3 country code (the panel unit).

year

Calendar year, 1990–2016.

contributions

Number of distinct UN peacekeeping operations the state contributes personnel to in that year (the count outcome).

democracy

Electoral-democracy index.

lgdppc

Log GDP per capita.

lpop

Log population.

milper

Military personnel, thousands.

majorpower

Major-power indicator (0/1).

region

World region.

Source

Contribution counts aggregated from the International Peace Institute Providing for Peacekeeping database (https://www.providingforpeacekeeping.org); covariates from V-Dem and the Correlates of War National Material Capabilities data. Counts and public covariates only; prepared for illustration.

Examples

data(peacekeeping)
table(peacekeeping$contributions == 0)

some <- subset(peacekeeping, iso3 %in% unique(iso3)[1:40])   # forty states keep the example short
## among these forty states every major power contributes, so majorpower would
## separate the participation equation (hurdle_cpb() warns when it does)
fit <- hurdle_cpb(contributions ~ lgdppc + milper, data = some,
                  participation = ~ democracy + lgdppc, fe = "iso3")
fit

Non-randomized PIT histogram for a fitted count model

Description

Draws the non-randomized probability integral transform histogram of Czado, Gneiting, and Held (2009). A well-calibrated model yields a flat histogram at height one (the reference line). A hump means the data are less dispersed than the fitted model predicts (underdispersed relative to it, as a Poisson fit to underdispersed counts shows), and a U shape means they are more dispersed.

Usage

pit_hist(
  fit,
  bins = 10,
  main = "PIT histogram",
  xlab = "PIT",
  ylab = "Relative frequency",
  ...
)

Arguments

Value

Invisibly, the vector of bin heights (normalized so that a calibrated model gives heights near one).

Examples

set.seed(1)
y <- rcpb(400, lambda = 3, alpha = 0.5)
fit <- cpb(y ~ 1, data = data.frame(y = y), truncated = FALSE, se = "none")
pit_hist(fit)

Predictions from a matched count-family fit

Description

Predictions from a matched count-family fit

Usage

## S3 method for class 'count_reg'
predict(
  object,
  newdata = NULL,
  type = c("response", "link", "prob"),
  at = NULL,
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Predictions from a CPB fit

Description

Predictions from a CPB fit

Usage

## S3 method for class 'cpb'
predict(
  object,
  newdata = NULL,
  type = c("response", "rate", "link", "ceiling", "prob"),
  at = NULL,
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Examples

set.seed(1); x <- rnorm(300)
N <- pmax(round(exp(1.6 + 0.5 * x) / 0.5), 1); y <- rbinom(300, N, 0.5)
fit <- cpb(y ~ x, data = data.frame(y = y, x = x)[y > 0, ], se = "none")
predict(fit, newdata = data.frame(x = 0), type = "ceiling")
predict(fit, newdata = data.frame(x = c(-1, 1)), type = "prob", at = 5)

Predictions from a fixed-effects CPB fit

Description

Predictions from a fixed-effects CPB fit

Usage

## S3 method for class 'cpb_fe'
predict(
  object,
  newdata = NULL,
  type = c("response", "rate", "link", "ceiling"),
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Predictions from a generalized event count fit

Description

Predictions from a generalized event count fit

Usage

## S3 method for class 'gec'
predict(
  object,
  newdata = NULL,
  type = c("response", "link", "prob"),
  at = NULL,
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Predictions from a fixed-effects GEC fit

Description

Predictions from a fixed-effects GEC fit

Usage

## S3 method for class 'gec_fe'
predict(
  object,
  newdata = NULL,
  type = c("response", "link"),
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Predict from a hurdle Poisson/NB fit

Description

Predict from a hurdle Poisson/NB fit

Usage

## S3 method for class 'hurdle_count'
predict(
  object,
  newdata = NULL,
  type = c("response", "participation", "intensity"),
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Predict from a hurdle-CPB fit

Description

Predict from a hurdle-CPB fit

Usage

## S3 method for class 'hurdle_cpb'
predict(
  object,
  newdata = NULL,
  type = c("response", "participation", "intensity"),
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Predict from a hurdle-GEC fit

Description

Predict from a hurdle-GEC fit

Usage

## S3 method for class 'hurdle_gec'
predict(
  object,
  newdata = NULL,
  type = c("response", "participation", "intensity"),
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Predict from a zero-inflated Poisson/NB fit

Description

Predict from a zero-inflated Poisson/NB fit

Usage

## S3 method for class 'zi_count'
predict(
  object,
  newdata = NULL,
  type = c("response", "count", "zero", "intensity"),
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Predict from a zi_cpb fit

Description

Predict from a zi_cpb fit

Usage

## S3 method for class 'zi_cpb'
predict(
  object,
  newdata = NULL,
  type = c("response", "zero", "intensity", "rate"),
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Predict from a zi_gec fit

Description

Predict from a zi_gec fit

Usage

## S3 method for class 'zi_gec'
predict(
  object,
  newdata = NULL,
  type = c("response", "zero", "intensity"),
  offset = NULL,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  ...
)

Arguments

Value

A numeric vector. With interval = "confidence", a matrix with columns fit, lwr and upr; with se.fit = TRUE, a list of fit (the vector or the matrix), se.fit and method (each row's interval source). Without se.fit, the source is the attribute method.

Simulate from the continuous parameter binomial

Description

Draws counts from the continuous parameter binomial (CPB) or its zero-truncated variant.

Usage

rcpb(n, lambda, alpha, truncated = FALSE, support = c("ceiling", "floor"))

Arguments

Details

King's CPB has the support ⁠0, ..., ceiling(lambda / (1 - alpha))⁠, which contains 1 for every positive rate. The floor support, ⁠0, ..., floor(lambda / (1 - alpha))⁠, is {0} for a rate whose ceiling is below 1, so its zero-truncated distribution does not exist. With truncated = TRUE such a draw is set to 1 and a warning reports how many there were: the stated parameters give that count probability zero, so data simulated this way are not data from the model (in a simulation study, keep every ceiling at 1 or above).

Value

An integer vector of counts.

Examples

set.seed(1)
table(rcpb(1000, lambda = 3, alpha = 0.5))

Residuals of a fitted model

Description

Response residuals (y minus the fitted mean) for every model class, and for the single-equation classes Pearson residuals, the response residual divided by the standard deviation of the fitted distribution (its exact variance; the conditional one of a zero-truncated fit).

Usage

## S3 method for class 'count_reg'
residuals(object, type = c("response", "pearson"), ...)

## S3 method for class 'gec'
residuals(object, type = c("response", "pearson"), ...)

## S3 method for class 'cpb_fe'
residuals(object, type = c("response", "pearson"), ...)

Arguments

Value

A numeric vector, in the rows of the data.

Simulate from the hurdle continuous parameter binomial

Description

Simulate from the hurdle continuous parameter binomial

Usage

rhurdle_cpb(n, lambda, alpha, p, support = c("ceiling", "floor"))

Arguments

Value

An integer vector: 0 for nonparticipants, a zero-truncated CPB draw otherwise.

Examples

set.seed(1)
y <- rhurdle_cpb(1000, lambda = 3, alpha = 0.5, p = 0.6)
mean(y == 0)

Hanging rootogram for a fitted count model

Description

Draws a Tukey hanging rootogram: bars for the observed frequencies hang from the curve of expected frequencies, both on the square-root scale. Bars that hang below the zero line mark counts the model under-predicts; bars that stop short mark counts it over-predicts. A fit with frequency weights counts each row as many times as its weight in both the observed and the expected frequencies.

Usage

rootogram(
  fit,
  kmax = NULL,
  main = "Rootogram",
  xlab = "Count",
  ylab = "sqrt(frequency)",
  ...
)

Arguments

Value

Invisibly, a data frame of count, observed, and expected frequencies.

Examples

set.seed(1)
y <- rcpb(400, lambda = 3, alpha = 0.5)
fit <- cpb(y ~ 1, data = data.frame(y = y), truncated = FALSE, se = "none")
rootogram(fit)

Simulate from the zero-inflated continuous parameter binomial

Description

Simulate from the zero-inflated continuous parameter binomial

Usage

rzicpb(n, lambda, alpha, pi, support = c("ceiling", "floor"))

Arguments

Value

An integer vector: a structural zero with probability pi, otherwise an (untruncated) CPB draw, which may itself be zero.

Examples

set.seed(1)
y <- rzicpb(1000, lambda = 3, alpha = 0.5, pi = 0.3)
mean(y == 0)

Proper scoring rules for a fitted count model

Description

Computes the mean logarithmic score and the ranked probability score (RPS) of a fitted model's predicted distribution against observed counts. Lower is better for both, and both are proper, so they are the natural way to compare the calibration of competing count models.

Usage

score(fit, newdata = NULL, kmax = NULL, weights = NULL, offset = NULL)

Arguments

Details

By default the scores are computed in sample (against the data the model was fit to, as means weighted by its frequency weights). Supply newdata to score a fitted model on held-out observations, or use cv_score() for a cross-validated score; an in-sample log score equals -logLik/n and does not penalize model complexity, so models of different size are better compared on the held-out or cross-validated score.

Value

A named numeric vector c(logscore, rps).

See Also

cv_score()

Examples

set.seed(1)
y <- rcpb(400, lambda = 3, alpha = 0.5)
fit <- cpb(y ~ 1, data = data.frame(y = y), truncated = FALSE, se = "none")
score(fit)

Simulate responses from a fitted underdisp model

Description

Draws nsim replicate response vectors from the fitted model, at the estimated parameters and the estimation data, in the format of stats::simulate(). Available for every model class in the package: cpb, cpb_fe, gec, gec_fe (through inheritance), count_reg, hurdle_count, zi_count, hurdle_cpb, hurdle_gec, zi_cpb, and zi_gec. Two-part models first draw the binary stage (participation or structural zero), then the count stage from its own distribution, so the replicates carry the model's full zero structure.

Usage

## S3 method for class 'cpb'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'cpb_fe'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'gec'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'count_reg'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'hurdle_cpb'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'hurdle_gec'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'hurdle_count'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'zi_cpb'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'zi_gec'
simulate(object, nsim = 1, seed = NULL, ...)

## S3 method for class 'zi_count'
simulate(object, nsim = 1, seed = NULL, ...)

Arguments

Details

The main consumer is simulated-residual diagnostics: DHARMa::createDHARMa(simulatedResponse = as.matrix(simulate(fit, 250)), observedResponse = y, fittedPredictedResponse = fitted(fit), integerResponse = TRUE) works for any fit in the family.

Value

A data frame with nsim integer columns, one row per observation, with a "seed" attribute.

Examples

set.seed(1); x <- rnorm(200)
N <- pmax(round(exp(1.4 + 0.4 * x) / 0.5), 1); y <- rbinom(200, N, 0.5)
fit <- cpb(y ~ x, data = data.frame(y = y, x = x)[y > 0, ], se = "none")
sims <- simulate(fit, nsim = 5)
colMeans(sims)

Summarize a CPB fit

Description

Summarize a CPB fit

Usage

## S3 method for class 'cpb'
summary(object, ...)

Arguments

Value

An object of class "summary.cpb" with the coefficient table, the dispersion parameter and its profile-likelihood interval (first-order, or calibrated when the fit went through calibrate_alpha()), the implied ceiling, fit statistics, and the likelihood-ratio test against a (zero-truncated) Poisson.

Tidy a fitted model (broom method)

Description

One row per coefficient with its estimate, standard error, z statistic and p-value, for every model class of the package. The two-part classes (⁠hurdle_*⁠, ⁠zi_*⁠) have a component column (the broom convention for multi-equation models) and the term prefixed by its component, so that modelsummary lays the table out without a shape argument. A term of a zero-inflation or participation equation that the data do not bound (separation; see zi_cpb()) has no standard error.

Usage

## S3 method for class 'cpb'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'cpb_fe'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'gec'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'count_reg'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'hurdle_cpb'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'hurdle_gec'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'hurdle_count'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'zi_cpb'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'zi_gec'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'zi_count'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

Arguments

Details

The scale column names what each coefficient multiplies: "log mean" for the mean-parameterized count families ("log untruncated mean" for a zero-truncated fit or a hurdle's intensity), the binary equation's link for a participation or zero-inflation row, and otherwise the family's own parameter ("log rate" for the CPB and GEC, "log lambda" for the rate-parameterized COM-Poisson, "log centrality" for the EWP, and so on; see irr.count). Coefficients on different scales are not comparable across the columns of a table; ame() puts every model's effect on the mean.

Value

A data frame with one row per coefficient.

Screen a count outcome for underdispersion

Description

Applies the diagnostic sequence developed in Bagozzi (2026): a marginal verdict from the conditional Pearson statistic and a through-origin score regression (which gives the direction of dispersion), a negative-binomial-versus-Poisson test for the overdispersion call, and—most importantly—an at-risk verdict on the positive counts (the units at risk of the event, in a hurdle's reading those that participate), benchmarked against a zero-truncated Poisson. The last step is what separates genuine underdispersion from the artifact of conditioning on Y>0. When the data are not zero-dominated it also fits the CPB, and it fits a generalized-Poisson soft-tail comparator. (The CPB's bound is tested on held-out counts with ceiling_exceedance(); in sample every feasible fit respects it.)

Usage

ud_screen(
  formula,
  data,
  run_cpb = TRUE,
  cpb_max_n = 3000,
  run_gp = TRUE,
  run_comp = TRUE,
  comp_max_par = 30,
  comp_max_n = 5000,
  ztp_threshold = c("calibrated", "bootstrap"),
  ztp_boot_B = 199L,
  cores = 1L,
  digits = 3,
  nb_boot = NULL,
  support = c("ceiling", "floor")
)

Arguments

Details

Runtime: the marginal and at-risk arms take seconds. ztp_threshold = "bootstrap" refits the zero-truncated Poisson ztp_boot_B times (use cores); the CPB comparator is fit only up to cpb_max_n rows, and the generalized-Poisson and COM-Poisson comparators only when the mean model carries at most comp_max_par parameters and at most comp_max_n rows, beyond which their rows print NA; each is a full model fit (the two soft-tail families have no concentrated fixed-effects path, so a dummy-heavy screen would take minutes).

Value

An object of class "ud_screen" with verdict_marginal, verdict_atrisk, the conditional and at-risk (ZTP-benchmarked) Pearson statistics, the NB-vs-Poisson LR test (p_nb_method records how its p-value was computed), a log-likelihood comparison, (when fit) the CPB alpha, and the over-conditioning guard state: atrisk_skipped (TRUE when the mean model nearly saturates the positive counts, so the at-risk statistic is not computed and the printout says why), sat_ratio (the fitted parameter share of the positives), and overconditioned (TRUE when that share reaches 0.10, the region where the calibrated threshold is anti-conservative; the printout then flags the verdict as diagnostic rather than probative and recommends ztp_threshold = "bootstrap").

References

King, G. (1989). Variance specification in event count models: From restrictive assumptions to a generalized estimator. AJPS 33(3):762-784.

Bagozzi, B. E. (2026). Revisiting Underdispersion in Political Science: Diagnostics and a Typology for Count Data. Companion manuscript.

See Also

cpb(), compare_dispersion()

Examples

set.seed(1); x <- rnorm(250)
N <- pmax(round(exp(1.4 + 0.4 * x) / 0.6), 1); y <- rbinom(250, N, 0.6)
ud_screen(y ~ x, data = data.frame(y = y, x = x))

Bootstrap covariance for a fixed-effects CPB fit

Description

The covariance of the covariate coefficients from the stored pairs/cluster bootstrap (se = "bootstrap" at fit time); NULL when no bootstrap was run.

Usage

## S3 method for class 'cpb_fe'
vcov(object, ...)

Arguments

Value

A covariance matrix, or NULL.

Zero-inflated count regression for the matched count families

Description

Fits a structural-zero mixture whose count component comes from any of the package's count families (see count_reg(), Families), the count analogue of zi_cpb(). The count component uses a log link (on the mean for the mean-parameterized families, on the rate \lambda for "compois") and the inflation probability a link set by link. Returned as a "zi_count" object that compare_models() and score() accept.

Usage

zi_count(
  formula,
  data,
  family = c("poisson", "negbin", "compois", "mpcmp", "genpois", "gammacount",
    "doublepois", "hyperpois", "ewp3", "ewp2"),
  zero = NULL,
  fe = NULL,
  zero_fe = NULL,
  link = c("logit", "probit", "cloglog"),
  offset = NULL,
  weights = NULL,
  se = c("analytic", "robust", "cluster", "none"),
  cluster = NULL
)

Arguments

Value

An object of class "zi_count"; ⁠$fitted.values⁠ is the marginal mean ⁠(1 - pi) E(Y | count component)⁠ that fitted() returns, and ⁠$mu⁠ the count component's natural parameter.

Separation

The fit checks its zero-inflation equation for separation as zi_cpb() does (see its Separation section), under its own link. It warns when the data put no useful bound on a covariate coefficient, and stores the check in the component separation. Such a coefficient's standard error is not meaningful either.

See Also

zi_cpb(), hurdle_count(), compare_models()

Examples

set.seed(2); n <- 300; x <- rnorm(n); z <- rnorm(n)
y <- ifelse(rbinom(n, 1, plogis(-0.5 + 0.8 * z)) == 1, 0L, rpois(n, exp(1 + 0.3 * x)))
zi_count(y ~ x, data.frame(y = y, x = x, z = z), family = "poisson",
         zero = ~ z, se = "none")

Zero-inflated continuous parameter binomial regression

Description

Fits the zero-inflated CPB. The model mixes a structural-zero process (a logistic model for the probability pi that a unit is a structural zero) with an untruncated CPB for the count, so a zero can arise either structurally or as a sampling zero from the CPB. Unlike hurdle_cpb(), the likelihood does not factorize, so the parameters are estimated by direct joint maximum likelihood. The mixture likelihood's optima differ with the structural share of zeros, so the search starts from points that span it: separate fits (a zero-truncated CPB on the positives and a logit of the zero indicator), the zero-inflated Poisson's coefficients, and a low-inflation start (a constant structural-zero probability of 0.02), followed by a BFGS pass; it is cross-checked against an expectation-maximization climb of the same observed-data likelihood, keeping the better optimum. The separation check (Separation) restarts the final polish from any better point it finds. An EM-only path (method = "em") is retained for comparison.

Usage

zi_cpb(
  formula,
  data,
  zero = NULL,
  fe = NULL,
  zero_fe = NULL,
  method = c("ml", "em"),
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  max.support = NULL,
  maxit = 200,
  tol = 1e-06,
  offset = NULL,
  weights = NULL,
  cores = 1L,
  support = c("ceiling", "floor")
)

Arguments

Details

The inflation equation uses a logit link. For a probit or cloglog inflation link, use zi_count() (Poisson/NB/COM-Poisson count), whose link argument covers the binary stage.

Value

An object of class "zi_cpb". Its component separation holds the separation check (see Separation).

Separation

A zero contributes ⁠log(pi + (1 - pi) f0)⁠ to the likelihood, which rises with pi, and a positive count contributes ⁠log(1 - pi) + log f(y)⁠, which falls with it. So a coefficient of the zero-inflation equation can grow without bound, and the likelihood with it, under three conditions. The zeros it pushes toward pi = 1 can go there freely, and so can the positives it pushes toward pi = 0. The count component's probability of a zero, f0, covers the zeros it pushes toward pi = 0. Unlike separation in a logistic regression, this needs no covariate that splits the zeros from the positives. So a hurdle_cpb() whose participation logit uses the same covariates can have a finite maximum while the mixture has none. The optimizer then stops somewhere along the divergent path. The zero equation's coefficients there are not estimates, and the log-likelihood approximates a supremum, so likelihood comparisons with the fit (likelihood ratio, AIC, Vuong) lack their usual reference distributions.

After the fit, the profile likelihood of each covariate coefficient in the zero equation is evaluated at the estimate. It is evaluated again at 10 and at 100 standard deviations of its covariate further from zero. Each time, the other zero-equation coefficients are re-maximized with the count component held at its estimate. The fit warns when a coefficient's profile 100 standard deviations out is less than qchisq(0.95, 1) / 2 below its value at the estimate. That coefficient's 95% profile interval then reaches past an effect of 100 on the log-odds per standard deviation, and the data put no useful bound on it. This covers a coefficient with no finite estimate and one whose finite maximum lies on a near-flat ridge, which a fit cannot tell apart. Holding the count component fixed can only lower the profile, so the full model's interval is at least as wide. Every profile can also reach the no-inflation boundary, where each structural-zero probability is zero and the zero equation drops out. A coefficient whose profile at the estimate lies within qchisq(0.95, 1) / 2 of that limit is unbounded because the data show little inflation, not by separation, so it is not flagged, and the fit is reported as near the boundary instead (zi_test() tests it). The intercept and the zero_fe dummies are not checked. Their limits are the no-inflation boundary (zi_test()) and the boundaries of single units (a unit with no positive count, whose structural-zero probability goes to one, as in a fixed-effects logit; with fe, such a unit is left out of the estimation, see fe). The component separation of the fit holds:

• diverging: the terms the data do not bound;

• profile: the profile steps;

• short_of_max: how far the best point the check found lies above the fit, with the count component unchanged. That point can come from the re-maximization at the estimate, when the fit is not a stationary point of the zero equation, or from a profile that rises above its value at the estimate, since the mixture likelihood is not concave and a profile can rise before it falls. A positive value means the optimizer stopped at least that far below the maximum;

• best: the zero-equation coefficients at that point (NULL when there is none);

• restarts, restart_gain: when the check finds a point at least 0.01 higher, the fit restarts its final polish there, with every parameter it estimates free, keeps the result only if it gains, and checks again, up to three times. These record the restarts kept and their total gain in log-likelihood. The other components describe the final fit;

• no_inflation, inflation_gain: whether the best point lies within qchisq(0.95, 1) / 2 of the no-inflation limit (the count component held), and its distance from it;

• n_extreme: the number of observations whose structural-zero probability is numerically 0 or 1, outside units that lie wholly at one edge.

zi_gec() and zi_count() run the same check.

Examples

set.seed(1)
n <- 800; x <- rnorm(n); z <- rnorm(n)
y <- rzicpb(n, lambda = exp(1.3 + 0.5 * x), alpha = 0.5, pi = plogis(-0.5 + 0.8 * z))
zi_cpb(y ~ x, data = data.frame(y = y, x = x, z = z), zero = ~ z)

Zero-inflated GEC (Katz-family) regression

Description

A structural-zero mixture with a GEC count component whose dispersion delta is estimated freely. Estimated by direct maximum likelihood from several starts, as in zi_cpb(): the single-equation GEC with a logit of the zero indicator for the inflation, at four dispersions; the zero-inflated Poisson's coefficients; and a low-inflation start. The best is polished, which avoids the degenerate pi -> 0 basin that coordinate-wise EM falls into.

Usage

zi_gec(
  formula,
  data,
  zero = NULL,
  zero_fe = NULL,
  se = c("none", "bootstrap"),
  B = 500,
  cluster = NULL,
  max.support = NULL,
  maxit = 200,
  tol = 1e-06,
  offset = NULL,
  weights = NULL,
  cores = 1L
)

Arguments

Value

An object of class "zi_gec".

Separation

The fit checks its zero-inflation equation for separation as zi_cpb() does (see its Separation section). It warns when the data put no useful bound on a covariate coefficient, and stores the check in the component separation.

See Also

zi_cpb(), gec(), hurdle_gec()

Examples

set.seed(4); n <- 300; x <- rnorm(n); z <- rnorm(n)
y <- ifelse(rbinom(n, 1, plogis(-0.5 + 0.8 * z)) == 1, 0L, rpois(n, exp(1 + 0.3 * x)))
zi_gec(y ~ x, data.frame(y = y, x = x, z = z), zero = ~ z)

Boundary-corrected test for zero-inflation (CPB vs ZI-CPB)

Description

A likelihood-ratio test of whether a count needs a structural-zero component. The (untruncated) CPB is nested in the zero-inflated CPB at a structural-zero probability of zero. Because that value lies on the boundary of the parameter space, the LR statistic follows a \tfrac12\chi^2_0 + \tfrac12\chi^2_1 mixture (Self and Liang 1987), which halves the naive \chi^2_1 p-value. The test uses an intercept-only structural-zero probability, so it is a clean single-parameter boundary test. Given two fitted models it refuses the cases where that reference distribution does not hold: a covariate-dependent structural-zero probability (several parameters on their boundary), a fixed-effects zi_cpb() (one count effect per unit, whose number grows with the panel, so the statistic has no fixed-dimension boundary reference), a separated zero equation, models fit to different data, a zero-truncated nest, and CPB models on different supports. It is not a test of a hurdle against a mixture, which are not nested. Compare those with compare_models().

Usage

zi_test(object, object2 = NULL, data = NULL, ...)

Arguments

Value

An object of class "zi_test" with the two log-likelihoods, the LR statistic, and the boundary-corrected p-value.

Examples

set.seed(1); x <- rnorm(500)
y <- rzicpb(500, lambda = exp(1.2 + 0.4 * x), alpha = 0.5, pi = 0.3)
zi_test(y ~ x, data = data.frame(y = y, x = x))