23 Poisson and negative binomial regression for count outcomes

23.1 Poisson regression

Poisson regression models the expected event count:

\[ Y_i\sim\operatorname{Poisson}(\mu_i), \]

\[ \log(\mu_i) = \beta_0+\beta_1X_{i1}+\cdots+\beta_qX_{iq}. \]

Exponentiating a coefficient gives a count ratio:

\[ \text{Count ratio}=e^\beta. \]

A ratio of 1 indicates no difference, 1.5 indicates a 50% higher expected count, and 0.80 indicates a 20% lower expected count.

If participants are observed for different lengths of time, the model should account for exposure time using an offset. Here, all participants have the same 12-week observation period, so no offset is required.

23.2 Fit the Poisson model

count_dat <- dat |>
  dplyr::select(ae_count_12w,treatment_arm,age_years,sex_at_birth,baseline_severity) |>
  drop_na() |>
  mutate(age10=(age_years-mean(age_years))/10)

count_formula <- ae_count_12w~treatment_arm+age10+sex_at_birth+baseline_severity

pois_fit <- glm(count_formula,family=poisson,data=count_dat)

23.3 Automated Poisson or negative binomial selection

Poisson regression assumes that the conditional variance is approximately equal to the conditional mean. When the variance is substantially larger, the data are overdispersed, and Poisson standard errors may be too small.

The Pearson dispersion statistic is used here to assess overdispersion.

pearson_chisq <- sum(residuals(pois_fit,type="pearson")^2)
df <- df.residual(pois_fit)
dispersion <- pearson_chisq/df
p_overdispersion <- pchisq(pearson_chisq,df,lower.tail=FALSE)

if(dispersion>1.5 && p_overdispersion<.05){
  count_fit <- MASS::glm.nb(count_formula,data=count_dat)
  count_method <- "Negative binomial regression"
} else {
  count_fit <- pois_fit
  count_method <- "Poisson regression"
}

tibble(
  dispersion=dispersion,
  p_overdispersion=p_overdispersion,
  selected_model=count_method
) |> kbl(digits=3,caption="Automated count-model selection")

A dispersion value near 1 is consistent with the Poisson model. In this teaching rule, a dispersion value above 1.5 together with \(p<0.05\) selects the negative binomial model. The value 1.5 is a practical teaching threshold, not a universal cutoff.

23.4 Adjusted count ratios

broom::tidy(count_fit,conf.int=TRUE,exponentiate=TRUE) |>
  filter(term!="(Intercept)") |>
  kbl(digits=3,caption="Adjusted count ratios")

A treatment count ratio of 1.6 means that the treatment group has an estimated 60% higher expected number of adverse events than Control over the same observation period, after adjustment for the other variables.

This should not be interpreted as 60% more participants experiencing an adverse event. The number of events and whether any event occurred are different outcomes.

23.5 Key takeaway

The outcome determines the model, while diagnostic checks help determine whether the chosen model needs refinement.