diagnose

diagnose()

diagnose() looks at a vector of counts and tells you three things that decide how to model it: whether it is overdispersed, whether it is bounded out of a fixed number of trials, and whether it has more zeros than its own dispersion explains. It then names the model family those signals imply. It does not fit the model for you; it hands you the diagnosis and points at the right shelf.

Usage

diagnose(x, k = NULL, group = NULL)

Argument	What it is
x	A non-negative integer vector, your count outcome. NAs are dropped.
k	The ceiling, if the count is "how many out of k" (items endorsed, questions correct). Leave it NULL for open-ended counts.
group	The design's condition factor, if the counts come from more than one group. Pass it so the dispersion test looks at variation within groups, not the spread created by a real group difference.

It returns a count_diagnosis object with a print method. The fields (vmr, prop_zero, bounded, dispersion_p_over, rec_code, and the rest) are all readable if you want to branch on them in code.

A bounded count: ACE scores out of 10

An adverse-childhood-experience score is how many of ten items a person endorses. Because it is bounded, the right home is the binomial family, not Poisson.

library(countkit)
set.seed(1)
ace <- sim_symptoms(n = 500, effect = 0, k = 10)$ace
diagnose(ace, k = 10)
#> Count diagnosis (n = 1000)
#>   mean = 1.53, var = 2.68, variance-to-mean ratio = 1.74
#>   zeros observed = 31.7% (Poisson predicts 21.5%; fitted NB predicts 31.8%)
#>   max = 8  (ceiling k = 10, bounded)
#>   dispersion ratio = 1.74 (overdispersion p < .001; underdispersion p = 1.000)
#>   >> Recommended: bounded-of-k -> binomial GLM with a dispersion check;
#>      beta-binomial or quasi-binomial under extra-binomial variation

Two things to read here. Because you passed k = 10, the tool sees the ceiling and routes to the binomial family. And look at the zero line: a third of the scores are zero, far past the 21.5% a Poisson would predict, but a fitted negative binomial expects almost exactly that many (31.8%). So the zeros are just what overdispersion produces, not a separate "never happened" process, and the tool does not falsely flag zero-inflation.

An open-ended, overdispersed count: verbal fluency

Fluency counts have no ceiling and are overdispersed. The tool sends them to the negative binomial.

set.seed(2)
fl <- sim_fluency(n = 300, effect = 0, mu = 15, vmr = 4)$y
diagnose(fl)
#> Count diagnosis (n = 600)
#>   mean = 14.99, var = 62.79, variance-to-mean ratio = 4.19
#>   zeros observed = 0.0% (Poisson predicts 0.0%; fitted NB predicts 0.1%)
#>   max = 46
#>   dispersion ratio = 4.18 (overdispersion p < .001; underdispersion p = 1.000)
#>   >> Recommended: overdispersed -> negative binomial (Poisson SEs will be too small)

Why group matters

This is the argument people forget, and it is the one that keeps the tool honest. If your counts come from two groups with genuinely different means, pooling them makes the outcome look overdispersed even when it is perfectly Poisson inside each group. Pass group and the dispersion test conditions on the group means.

set.seed(4)
d <- sim_fluency(n = 500, effect = log(2), mu = 15, vmr = 1)  # Poisson within group

diagnose(d$y)                    # pooled: fooled
#>   variance-to-mean ratio = 3.57
#>   dispersion ratio = 3.57 (overdispersion p < .001; underdispersion p = 1.000)
#>   >> Recommended: overdispersed -> negative binomial (Poisson SEs will be too small)

diagnose(d$y, group = d$group)   # conditional: not fooled
#>   dispersion ratio = 0.99 (overdispersion p = .599; underdispersion p = .401)
#>   >> Recommended: ~equidispersed, high mean -> Poisson fine; Gaussian approx. often acceptable

Same data, opposite advice. The pooled call sees a variance-to-mean ratio of 3.57 and says negative binomial; the conditional call sees that, within groups, the ratio is 0.99 and correctly says the counts are Poisson. Whenever you have a design, pass the design.

What it can recommend

The rec_code field carries a short machine-readable tag, so you can branch on the recommendation in code:

rec_code	When	Points to
betabinom	Bounded out of k, with extra-binomial variance	Beta-binomial or quasi-binomial
zi_binomial	Bounded and genuinely zero-inflated	Zero-inflated binomial or beta-binomial hurdle
nb	Open-ended and overdispersed	Negative binomial
qpois	Underdispersed (variance below the mean)	Quasi-Poisson or Conway-Maxwell-Poisson
hurdle	More zeros than a negative binomial predicts	Hurdle or zero-inflated model
poisson	Equidispersed	Poisson (a Gaussian approximation is often fine at high means)

The tool stops at the recommendation. Fitting the model it points to, and reporting coverage rather than only a p value, is the next step, and it is one line in each case (for example MASS::glm.nb(y ~ group) for the negative binomial).