Why Compare Count Models?
Repeated clinical counts can show several features at once: average event rates may change with treatment or exposure, observations from one person may be correlated, and the data may contain more variability or more zeros than a Poisson model predicts. A model comparison is useful only if it asks which of these features the candidate models represent, and whether their differences improve the description and prediction of the observed data.
This tutorial compares three count-response choices in an NLME setting: Poisson, negative binomial, and zero-inflated. Each can include a subject-specific random effect and an exposure offset. The comparison should use the likelihood for the observed repeated data, not merely the fit to counts conditional on estimated random effects. Predicted count distributions and rate effects provide important context for interpreting a numerical model-selection result.
Candidate Distributions and Their Assumptions
Under a Poisson model, \(Y_{ij}\mid b_i\) has mean and variance \(\mu_{ij}\). The exposure term enters multiplicatively, so \(\log(E_{ij})\) is an offset with coefficient fixed at one. This represents a rate that is constant within the interval after accounting for the modeled covariates and random effect.
A common negative-binomial (NB2) model retains mean \(\mu_{ij}\), but adds a dispersion parameter \(\kappa>0\). Its conditional variance is \(\mu_{ij}+\mu_{ij}^{2}/\kappa\), which exceeds the Poisson variance for finite \(\kappa\). The negative binomial can therefore represent extra count variation beyond that explained by the random effects and covariates.
A zero-inflated model assumes two possible sources of an observed zero: a structural-zero state with probability \(\pi_{ij}\), or a count-generating state with probability \(1-\pi_{ij}\). In the latter state, the count follows a Poisson or negative-binomial distribution. For a Poisson count component with mean \(\mu_{ij}\), the zero probability is \(\pi_{ij}+(1-\pi_{ij})\exp(-\mu_{ij})\), while the probability of a positive count \(y\) is \((1-\pi_{ij})f(y\mid\mu_{ij})\). The inflation probability may have its own covariate model. A structural-zero interpretation should be scientifically plausible; a high observed zero fraction alone does not establish a distinct structural-zero process.
In all three cases, a subject-level likelihood contribution integrates over the random effect. If \(f(y_{ij}\mid b_i)\) is the appropriate conditional probability mass function, then the contribution for subject \(i\) is:
Here \(\phi\) is the normal density of the random effect. In general, the integral of the product is not the product of probabilities evaluated at \(b_i=0\), nor is it the product evaluated at an estimated individual random effect. The marginal likelihood accounts for uncertainty in the random effects and is the appropriate basis for comparing these fitted NLME models.
Worked Example: Comparing Predicted Probabilities at a Common Mean
Consider one interval with exposure one and a conditional mean of 1.5 events. Compare a Poisson model, an NB2 model with \(\kappa=1.5\), and a zero-inflated Poisson model with structural-zero probability \(\pi=0.25\). These are illustrative conditional distributions at a fixed random effect, not population-averaged predictions.
For the Poisson model, \(P(Y=0)=e^{-1.5}=0.2231\). For \(Y=2\), the probability is \(e^{-1.5}(1.5^2/2!)=0.2510\).
For NB2, the success probability in the negative-binomial parameterization is \(p=\kappa/(\kappa+\mu)=1.5/(1.5+1.5)=0.5\). Thus \(P(Y=0)=p^\kappa=0.5^{1.5}=0.3536\). For \(Y=2\):
For the zero-inflated Poisson model, the count-generating state must have mean \(1.5/(1-0.25)=2\) to keep the overall mean at 1.5. Therefore \(P(Y=0)=0.25+0.75e^{-2}=0.3515\), and \(P(Y=2)=0.75e^{-2}(2^2/2!)=0.2030\).
The models have the same overall mean but different predicted probabilities for zero and for two events. The NB2 and zero-inflated models place substantially more probability on zero than the Poisson model here, but their distributions are not interchangeable: the NB2 model represents extra count variation, whereas the zero-inflated model includes a separate zero-generating state.
Using Likelihood and Information Criteria
For candidate models fitted to the same response observations, exposure definitions, and analysis set, a higher maximized log likelihood indicates a closer likelihood fit, but a more flexible model will generally fit at least as well. Akaike’s information criterion (AIC) balances fit against the number of estimated parameters: \(\mathrm{AIC}=-2\ell+2k\), where \(\ell\) is the maximized marginal log likelihood and \(k\) is the number of estimated parameters. Lower AIC is preferred among the candidates considered, but it is not a test that proves one model is true.
Worked Example: Comparing Three Hypothetical NLME Fits
Suppose a hypothetical clinical dataset contains repeated event counts from 24 subjects over four intervals, with differing interval exposures. Three models use the same log-rate covariates, subject random-intercept structure, and analysis observations; they differ only in count distribution. Maximum likelihood fitting gives the following illustrative results.
| Model | Marginal log likelihood \(\ell\) | Parameters \(k\) | AIC |
|---|---|---|---|
| Poisson | -42.8 | 5 | 95.6 |
| Negative binomial | -36.1 | 6 | 84.2 |
| Zero-inflated Poisson | -35.4 | 7 | 84.8 |
For the negative-binomial model, \(\mathrm{AIC}=-2(-36.1)+2(6)=72.2+12=84.2\). For the zero-inflated model, \(-2(-35.4)+2(7)=70.8+14=84.8\); for Poisson, \(-2(-42.8)+2(5)=85.6+10=95.6\). The negative-binomial model has the lowest AIC, but it is only 0.6 AIC units below the zero-inflated model. That small difference does not support a decisive distinction between those two candidates. Both fit substantially better by AIC than Poisson in this illustrative comparison.
The log likelihood improvement from Poisson to negative binomial is \(2[-36.1-(-42.8)]=13.4\) on the likelihood-ratio scale. It may be tempting to use a routine one-degree-of-freedom likelihood-ratio test, but the Poisson limit corresponds to a boundary value of the NB dispersion parameter. Standard regular chi-square reference assumptions may therefore fail. AIC is used here as a comparative criterion, not as a significance test. The zero-inflated model is also not nested in the negative-binomial model in the usual way, so the difference between them is not interpreted through that likelihood-ratio calculation.
Short Justification: Why the Random-Effects Integral Matters
Conditional independence does not imply marginal independence. For two observations from the same subject, their joint probability is obtained by averaging the conditional product over the shared random effect:
The same \(b_i\) appears in both factors. If instead one multiplied the two separately averaged probabilities, the calculation would be \(\{\int f(y_{i1}\mid b)\phi(b)\,db\}\{\int f(y_{i2}\mid b)\phi(b)\,db\}\), which corresponds to drawing an independent random effect for each interval. That is a different model and removes the within-subject dependence induced by the shared effect. This is why model-comparison likelihoods must use the specified subject-level integration.
Worked Example: A Repeated-Count Likelihood Contribution
For an arithmetic illustration, suppose two counts from one hypothetical subject are 0 and 2, and the fitted conditional Poisson means at a specified random-effect value are 1.5 and 4.5. Conditional independence gives the product of the two Poisson probabilities:
The arithmetic uses \(4.5^2/2=10.125\) and \(e^{-6}\approx0.00247875\), giving \(0.0250974\). This is a conditional probability product at the specified random-effect value. For the actual NLME subject likelihood, the product must be integrated over the distribution of \(b_i\); the displayed value alone is not the marginal likelihood contribution.
From Model Coefficients to Clinical Rates
With a log link and exposure offset, a treatment coefficient \(\beta_T\) compares rates at the same exposure and other covariate values. Its exponent, \(\exp(\beta_T)\), is a rate ratio. A positive coefficient gives a ratio above one; a negative coefficient gives a ratio below one. The random intercept changes a subject’s rate multiplicatively, and the population prediction averages over its distribution rather than setting it to zero.
Worked Example: Interpreting a Treatment Rate Ratio
Suppose an NB2 NLME fit to hypothetical recurrent symptom counts uses exposure in patient-days and estimates the treatment coefficient as \(\hat{\beta}_T=-0.2877\), with a 95% confidence interval from \(-0.5108\) to \(-0.06188\). The estimated rate ratio is \(\exp(-0.2877)=0.7500\); exponentiating the interval endpoints gives \((0.6000,0.9400)\).
At the same exposure and covariate values, the estimated symptom-event rate in the treatment group is 0.75 times the comparison-group rate, a 25% lower estimated rate. The interval corresponds to rate ratios from 0.60 to 0.94, or rates 6% to 40% lower. Under the model and the usual repeated-sampling interpretation, this 95% confidence procedure has 95% coverage in repeated samples; it does not mean there is a 95% probability that the fixed rate ratio lies in this particular interval. The contrast is an adjusted rate association; a causal interpretation requires an appropriate study design and assumptions beyond the count model.
Diagnostics, Alternatives, and Limitations
AIC should be read alongside model diagnostics. Compare observed and predicted count frequencies, including the proportions with zero, one, and larger counts. Examine predictions by exposure, time, treatment, and subject-level pattern. For a fitted NLME model, predictive checks should reflect both the count distribution and the random-effects structure; otherwise the simulated repeated-count patterns will not be comparable to the data.
- Poisson inadequacy: residual or predictive patterns may show more variability or more extreme counts than predicted. A random intercept can already create marginal overdispersion, so diagnose the full fitted model rather than comparing the raw variance with the raw mean alone.
- Zero-inflation concerns: the observed zero fraction may exceed predictions, or the model may predict too many zeros in some exposure groups. A high zero fraction does not by itself prove a structural-zero subgroup; check whether covariates, exposure, heterogeneity, or a different count distribution can explain it.
- Random-effects and optimization concerns: inspect convergence, parameter uncertainty, and sensitivity to starting values. A dispersion estimate near the Poisson limit or an inflation probability near a boundary can make model distinctions weak and numerical estimation unstable.
- Exposure and mean structure: check that the offset is appropriate and that relevant time trends or covariates have not been omitted. A distributional change cannot reliably repair a wrongly specified rate trend.
A hurdle model is a relevant alternative when the process determining any zero versus a positive count is scientifically distinct from the process determining the positive count distribution. It is not simply another name for zero inflation: in a hurdle model, the positive-count component is truncated at zero. The choice among Poisson, NB, zero-inflated, and hurdle models therefore concerns both statistical fit and the meaning of the data-generating assumptions.
Model comparison is limited to the candidates, covariates, and random-effects structures actually fitted. Similar AIC values indicate that the data do not clearly separate those candidates by this criterion; they do not establish that the models are equivalent. Predictions can also be sensitive to whether they are conditional on a typical subject (\(b_i=0\)) or averaged over the subject population. State which prediction is being reported, and retain the exposure units when describing a rate.
Check Your Understanding
Answer each question, then compare with the solutions below.
- Why is a subject’s marginal likelihood contribution not generally the product of probabilities evaluated at \(b_i=0\)?
- Two models have AIC values 84.2 and 84.8. What is the appropriate interpretation of the difference?
- An NB2 model has mean 3 and \(\kappa=2\). What is its conditional variance?
- A treatment coefficient is \(\log(0.80)\) in a log-rate model. Interpret the rate ratio.
Solutions. (1) The likelihood contribution integrates the conditional product over the shared random-effect distribution; evaluation at zero ignores random-effect uncertainty and heterogeneity. (2) The first model has the lower AIC, but a difference of 0.6 is small and does not decisively distinguish the candidates. (3) The variance is \(3+3^2/2=7.5\) counts squared for the modeled interval. (4) The estimated treatment rate is 0.80 times the comparison rate at equal exposure and covariates, or 20% lower; this is an adjusted rate contrast, not by itself a causal conclusion.