Negative-Binomial Models for Overdispersed Counts in Joint Models

Why Use a Negative-Binomial Count Model?

Repeated counts often vary more than a Poisson model allows. For a Poisson random variable with mean $\mu$, the variance is also $\mu$. In longitudinal data, however, counts with similar fitted means can differ substantially in their variability. If that extra variation is ignored, standard errors and uncertainty statements can be unreliable, and predictions may be too concentrated.

A negative-binomial longitudinal submodel offers one way to represent this extra conditional variation. It does not, by itself, account for every source of heterogeneity between people: subject-specific random effects can represent persistent differences, while the negative-binomial dispersion parameter represents additional variation remaining after conditioning on those effects and measured covariates. Keeping these sources distinct is central to specifying and interpreting the model.

Definition: In the NB2 parameterization, a count $Y$ has conditional mean $\mu>0$ and dispersion parameter $\alpha\geq 0$, with $$ \operatorname{Var}(Y\mid \mu)=\mu+\alpha\mu^2. $$ When $\alpha>0$, this variance exceeds the Poisson variance $\mu$. The Poisson model is the limiting case as $\alpha$ approaches zero.

Longitudinal NB2 Submodel in a Joint Model

Let $Y_{ij}$ be the count for subject $i$ at measurement occasion $j$, and let $e_{ij}>0$ be the exposure during that occasion, such as person-time at risk. A common model uses a log link:

Y_{ij}\mid b_i \sim \operatorname{NB2}(\mu_{ij},\alpha),\qquad \log(\mu_{ij})=\log(e_{ij})+x_{ij}^{\mathsf T}\beta+z_{ij}^{\mathsf T}b_i.

Here, $x_{ij}$ is a vector of observed predictors, $\beta$ is a vector of fixed-effect coefficients, $z_{ij}$ is a vector of random-effect design values, and $b_i$ is a subject-specific random-effects vector, commonly modeled as $b_i\sim N(0,D)$. The offset $\log(e_{ij})$ makes the expected count proportional to exposure when the other predictors are held fixed. Without an exposure offset, the linear predictor describes a mean count for the observation interval represented by the model, not automatically a rate.

Conditional on $b_i$ and the covariates, the count has mean $\mu_{ij}$ and variance $\mu_{ij}+\alpha\mu_{ij}^2$. A typical specification also assumes conditional independence of the repeated count observations given the random effects and covariates. In a joint model, the event-time component can share some or all of $b_i$ with the longitudinal component. The shared effects allow the two processes to be associated through subject-level latent characteristics; the NB2 dispersion does not replace that linkage. The joint model also needs its own assumptions about the event process and how it depends on shared effects.

The random effects and NB2 dispersion address different aspects of variability. Conditional on $b_i$, $\alpha$ controls extra count variation beyond the Poisson variance. Across subjects, variation in $b_i$ also changes expected counts. Consequently, the unconditional variance across subjects can exceed the conditional NB2 variance. It is not generally correct to interpret $\alpha$ as a summary of all between- and within-subject heterogeneity.

Why the NB2 Variance Has a Quadratic Term

The NB2 distribution can be represented as a Poisson count whose rate itself varies. In this representation, let $\Lambda$ be a gamma-distributed rate with mean $\mu$ and variance $\alpha\mu^2$, and let $Y\mid\Lambda\sim\operatorname{Poisson}(\Lambda)$. The law of total variance gives a direct justification for the NB2 variance:

\begin{aligned} \operatorname{Var}(Y) &=E\{\operatorname{Var}(Y\mid\Lambda)\} +\operatorname{Var}\{E(Y\mid\Lambda)\}\\ &=E(\Lambda)+\operatorname{Var}(\Lambda)\\ &=\mu+\alpha\mu^2. \end{aligned}

The first term is the usual Poisson variation conditional on the latent rate. The second is variation in that rate. For the gamma distribution in this representation, the shape is $1/\alpha$ and the scale is $\alpha\mu$; their product is the mean $\mu$, and the shape times the squared scale is the variance $\alpha\mu^2$. This also explains why the NB2 variance grows faster than the mean as $\mu$ increases.

Worked Example: Comparing Poisson and NB2 Variances

Suppose a hypothetical measurement occasion has fitted conditional mean $\mu=4$ counts and NB2 dispersion $\alpha=0.25$. Compare the conditional variances and calculate the probability of a zero count under NB2.

Under a Poisson model, the variance equals the mean, so the variance is $4$ counts squared. Under NB2,

\operatorname{Var}(Y\mid b_i)=4+0.25(4^2)=4+4=8.

Thus, for this same conditional mean, the NB2 variance is $8$ counts squared, twice the Poisson variance. Under this parameterization, the NB2 size is $1/\alpha=4$. Its zero-count probability is

P(Y=0\mid b_i)=(1+\alpha\mu)^{-1/\alpha} =(1+0.25(4))^{-4}=2^{-4}=0.0625.

This probability is 6.25% for the stated conditional mean and dispersion. The example illustrates the model’s distributional implication; it does not establish that NB2 is suitable for any particular data set.

Interpreting Rates and Expected Counts

With a log link and an exposure offset, the linear predictor apart from $\log(e_{ij})$ is the log rate. Exponentiating it gives the conditional rate, in expected counts per unit of exposure. A coefficient for a one-unit increase in a predictor is a log rate ratio, so exponentiating the coefficient gives the multiplicative rate ratio, holding the other predictors and random effects fixed. To obtain an expected count over a specified exposure period, multiply the rate by that exposure.

Worked Example: A Treatment Rate Ratio and Expected Count

Consider a hypothetical model with a baseline rate of $2$ counts per person-month at time zero. Let treatment be coded $1$ for treated and $0$ otherwise, with treatment coefficient $-0.30$, and let time $t$, measured in years, have coefficient $0.10$. For illustration, take the subject-specific random effect to be zero. The model is

\log(\mu_{ij})=\log(e_{ij})+\log(2)-0.30\,\text{treatment}_{ij}+0.10t_{ij}.

The treatment rate ratio at the same time and exposure is

\exp(-0.30)=0.7408.

Thus, the model predicts a rate about $0.741$ times as high for treated observations as for untreated observations, conditional on the stated predictors and random effect. This is a rate comparison, not a comparison of the NB2 variances.

At year $2$, the untreated rate is $2\exp(0.10(2))=2\exp(0.20)\approx2.443$ counts per person-month. The treated rate is $2\exp(-0.30+0.20)=2\exp(-0.10)\approx1.810$ counts per person-month. For a treated observation with exposure $1.5$ person-months, the expected count is

1.5\{2\exp(-0.10)\}\approx1.5(1.8097)=2.715\text{ counts}.

The value $2.715$ is an expected count, not a prediction that the observed count must be close to a non-integer value. The count outcome remains integer-valued and variable.

Estimating and Checking Overdispersion

In an actual longitudinal joint model, $\alpha$ is estimated together with the fixed effects, random-effects distribution, and event-process parameters under the specified joint likelihood. A simple method-of-moments calculation can build intuition, but it does not account for covariates, random effects, or dependence between repeated observations, and should not be substituted for fitting the joint model.

Worked Example: A Simple Moment-Based Dispersion Calculation

Suppose a hypothetical set of five independent counts from a common NB2 distribution is $0,1,2,3,6$. This deliberately simplified example has no repeated-subject structure or covariates. The sample mean is $\bar{x}=12/5=2.4$. The sum of squared observations is $50$, so the sample variance is

s^2=\frac{50-12^2/5}{5-1} =\frac{50-28.8}{4}=5.3.

Equating the NB2 variance to the sample variance and replacing the population mean by $\bar{x}$ gives the moment estimate

\hat{\alpha}_{\mathrm{mom}} =\frac{s^2-\bar{x}}{\bar{x}^2} =\frac{5.3-2.4}{2.4^2} =\frac{2.9}{5.76}\approx0.5035.

The estimate is positive because the sample variance exceeds the sample mean. With only five observations it is highly uncertain. In longitudinal data, the same calculation on pooled raw counts could mistake subject heterogeneity or changing fitted means for NB2 dispersion; the joint model’s covariates and random effects must be considered.

Useful diagnostics ask whether the fitted distribution reproduces variability at the level being modeled. Compare observed and fitted count distributions across relevant ranges of fitted means, and examine residual patterns over time and by exposure or covariate groups. Large residual spread relative to what the fitted model predicts can indicate inadequate dispersion modeling, a misspecified mean trajectory, omitted predictors, or dependence not captured by the random effects. A dispersion estimate near zero suggests little evidence for NB2 variation beyond the Poisson component, though uncertainty and model fit still matter.

A negative-binomial model does not automatically solve every mismatch in the count distribution. In particular, poor fit among zero counts may signal that the overall model or mean structure is inadequate; specialized models for excess zeros are a distinct alternative. Also check whether exposure is correctly measured and whether proportionality to exposure is plausible. If a two-unit exposure interval does not reasonably have twice the expected count of a one-unit interval, a simple log-exposure offset may be inappropriate.

Choosing Between Poisson and Negative Binomial

Feature	Poisson longitudinal submodel	NB2 longitudinal submodel
Conditional variance	$\mu_{ij}$	$\mu_{ij}+\alpha\mu_{ij}^2$
Extra conditional dispersion	Not represented beyond the Poisson variance	Represented by $\alpha\geq0$
Subject heterogeneity	Can be modeled using shared or subject-specific random effects	Can also be modeled using random effects; $\alpha$ is additional to them
Interpretation of log-link coefficients	Conditional log rate ratios	Conditional log rate ratios, with the same rate interpretation

Use a Poisson submodel when its conditional variance structure is credible after accounting for predictors and random effects. NB2 is a reasonable extension when meaningful extra conditional variation remains and its quadratic mean-variance relationship fits the data. The choice should not be made solely because the observed raw variance exceeds the raw mean: changing means over time and variation between subjects can produce that pattern even when the conditional count distribution is Poisson.

The NB2 model is also not the only negative-binomial parameterization. For example, an NB1 formulation uses a variance that grows linearly with the mean rather than quadratically. The parameter and variance interpretation must therefore be stated explicitly; a reported dispersion value is not interpretable without knowing the chosen parameterization.

Common Mistakes and Limitations

Confusing dispersion with the mean. $\alpha$ does not add a fixed amount to the expected count. It changes the variance through $\alpha\mu^2$.
Calling every excess variance “NB2 dispersion.” First account for covariates, time trends, exposure, and subject-level random effects. These can all create marginal variability.
Interpreting a rate as a count. A rate has units of counts per exposure unit; an expected count requires multiplying by the specified exposure.
Assuming a coefficient is a causal effect. An exponentiated treatment coefficient is a model-based conditional rate ratio. Causal interpretation requires an appropriate study design and additional assumptions.
Assuming NB2 guarantees a good fit. Incorrect mean structure, unusual zero frequencies, or unmodeled dependence can remain after adding $\alpha$.

The model’s predictions are conditional on the specified covariates and random effects. A population-average prediction instead averages over the random-effects distribution, and in a nonlinear log-link model that average is generally not obtained simply by setting $b_i=0$. Such predictions require integrating over the random-effects distribution. In a joint model, predictions can also depend on event-process information and the shared-effects structure.

Key takeaway: An NB2 longitudinal submodel retains the log-rate interpretation of a Poisson model but replaces the conditional variance $\mu$ with $\mu+\alpha\mu^2$. Interpret $\alpha$ as extra conditional dispersion, not as a substitute for subject random effects, and convert predicted rates to expected counts using the relevant exposure.

Check Your Understanding

Try each question, then check the solutions below.

If $\mu=5$ and $\alpha=0.20$, what is the NB2 conditional variance, and how does it compare with the Poisson variance?
A log-rate coefficient is $0.40$ for a one-unit predictor increase. What is the corresponding rate ratio, approximately?
A fitted rate is $1.492$ counts per person-week. What expected count does the model imply over three person-weeks?
Why can a pooled sample variance larger than its pooled mean fail to demonstrate NB2 dispersion in longitudinal data?
What does the offset $\log(e_{ij})$ assume about the relationship between exposure and expected count?

Solutions. (1) The variance is $5+0.20(5^2)=10$ counts squared, compared with Poisson variance $5$. (2) The rate ratio is $\exp(0.40)\approx1.492$, or about 1.49 times the rate per one-unit increase, conditional on the other model terms. (3) The expected count is $3(1.492)=4.476$ counts over three person-weeks. (4) Pooled variation can reflect changing means, subject random effects, or other heterogeneity, not only conditional NB2 dispersion. (5) It assumes expected count is proportional to exposure, holding the remaining predictors and random effects fixed.

Feature	Poisson longitudinal submodel	NB2 longitudinal submodel
Conditional variance	\(\mu_{ij}\)	\(\mu_{ij}+\alpha\mu_{ij}^2\)
Extra conditional dispersion	Not represented beyond the Poisson variance	Represented by \(\alpha\geq0\)
Subject heterogeneity	Can be modeled using shared or subject-specific random effects	Can also be modeled using random effects; \(\alpha\) is additional to them
Interpretation of log-link coefficients	Conditional log rate ratios	Conditional log rate ratios, with the same rate interpretation