Models for Count Outcomes

Poisson regression and variations

Acknowledgements

This content is adapted from:

• Dobson and Barnett (2018), Chapter 9
• Vittinghoff et al. (2012), Chapter 8

Configuring R

Functions from these packages will be used throughout this document:

library(conflicted) # check for conflicting function definitions
# library(printr) # inserts help-file output into markdown output
library(rmarkdown) # Convert R Markdown documents into a variety of formats.
library(pander) # format tables for markdown
library(ggplot2) # graphics
library(ggfortify) # help with graphics
library(dplyr) # manipulate data
library(tibble) # `tibble`s extend `data.frame`s
library(magrittr) # `%>%` and other additional piping tools
library(haven) # import Stata files
library(knitr) # format R output for markdown
library(tidyr) # Tools to help to create tidy data
library(plotly) # interactive graphics
library(dobson) # datasets from Dobson and Barnett 2018
library(parameters) # format model output tables for markdown
library(haven) # import Stata files
library(latex2exp) # use LaTeX in R code (for figures and tables)
library(fs) # filesystem path manipulations
library(survival) # survival analysis
library(survminer) # survival analysis graphics
library(KMsurv) # datasets from Klein and Moeschberger
library(parameters) # format model output tables for
library(webshot2) # convert interactive content to static for pdf
library(forcats) # functions for categorical variables ("factors")
library(stringr) # functions for dealing with strings
library(lubridate) # functions for dealing with dates and times
library(broom) # Summarizes key information about statistical objects in tidy tibbles
library(broom.helpers) # Provides suite of functions to work with regression model 'broom::tidy()' tibbles

Here are some R settings I use in this document:

rm(list = ls()) # delete any data that's already loaded into R

conflicts_prefer(dplyr::filter)
ggplot2::theme_set(
  ggplot2::theme_bw() + 
        # ggplot2::labs(col = "") +
    ggplot2::theme(
      legend.position = "bottom",
      text = ggplot2::element_text(size = 12, family = "serif")))

knitr::opts_chunk$set(message = FALSE)
options('digits' = 6)

panderOptions("big.mark", ",")
pander::panderOptions("table.emphasize.rownames", FALSE)
pander::panderOptions("table.split.table", Inf)
conflicts_prefer(dplyr::filter) # use the `filter()` function from dplyr() by default
legend_text_size = 9
run_graphs = TRUE

1 Introduction

This chapter presents models for count data outcomes. With covariates, the event rate \(\lambda\) becomes a function of the covariates \(\tilde{X}= (X_1, \dots,X_n)\). Typically, count data models use a \(\operatorname{log}\mathopen{}\left\{\right\}\mathclose{}\) link function, and thus an \(\operatorname{exp}\mathopen{}\left\{\right\}\mathclose{}\) inverse-link function. That is:

\[ \begin{aligned} \operatorname{E}[Y | \tilde{X}= \tilde{x}, T = t] &= \mu(\tilde{x},t) \\ \mu(\tilde{x},t) &= \lambda(\tilde{x})\cdot t \\ \lambda(\tilde{x}) &= \operatorname{exp}\mathopen{}\left\{\eta(\tilde{x})\right\}\mathclose{} \\ \eta(\tilde{x}) &= \tilde{x}'\tilde \beta = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p \end{aligned} \tag{1}\]

\(T=t\) is called the exposure magnitude and has a special role in this model.

Exercise 1 Where have we seen a relationship like \[\mu = \lambda \cdot t\] before?

Solution 1. The relationship \[\mu = \lambda \cdot t\] in count regression models is analogous to the relationship \[\mu = n \pi\] in Binomial models.

We can also think of \(t\) as a special part of the linear component:

\[ \begin{aligned} \operatorname{log}\mathopen{}\left\{\operatorname{E}[Y | \tilde{X}= \tilde{x},T=t]\right\}\mathclose{} &= \operatorname{log}\mathopen{}\left\{\mu(\tilde{x})\right\}\mathclose{}\\ &=\operatorname{log}\mathopen{}\left\{\lambda(\tilde{x}) \cdot t\right\}\mathclose{}\\ &=\operatorname{log}\mathopen{}\left\{\lambda(\tilde{x})\right\}\mathclose{} + \log{t}\\ &=\operatorname{log}\mathopen{}\left\{\operatorname{exp}\mathopen{}\left\{\eta(\tilde{x})\right\}\mathclose{}\right\}\mathclose{} + \log{t}\\ &=\eta(\tilde{x}) + \log{t}\\ &=\tilde{x}'\tilde\beta + \log{t}\\ &=(\beta_0 +\beta_1 x_1+\dots + \beta_p x_p) + \log{t}\\ \end{aligned} \]

In contrast with the other covariates (represented by \(\tilde{X}\)), \(t\) enters this expression with a \(\log{}\) transformation and without a corresponding \(\beta\) coefficient; in other words, \(\operatorname{log}\mathopen{}\left\{t\right\}\mathclose{}\) is an offset term.

Exercise 2 What are the units of \(\mu\) in Equation 1?

Solution 2. \(\mu\) is the mean of \(Y\), and \(Y\) is a count, so \(\mu\) is also a count; for example:

• 3.1 cyclones,
• 10.23 ER visits
• 15.01 infections

Exercise 3 What are the units of \(\lambda\) in Equation 1?

Solution 3. \(\lambda = \mu/t\), so \(\lambda\) is a rate of counts per unit of \(t\). For example:

• 3.1 cyclones per year
• 2.023 ER visits per 10 person-years
• 15.01 infections per 1000 person-years at risk

2 Interpreting Poisson regression models

Applying the general procedure for interpreting a regression coefficient to the log-rate linear predictor \(\eta(\tilde{x}) = \operatorname{log}\mathopen{}\left\{\lambda(\tilde{x})\right\}\mathclose{}\) shows that \(\beta_j\) is the difference in log-rates per one-unit increase in \(x_j\), holding all other covariates fixed (and equal to their reference values when interactions are present).

Differences on the log-rate scale become ratios on the rate scale, because

\[\operatorname{exp}\mathopen{}\left\{a-b\right\}\mathclose{} = \frac{\operatorname{exp}\mathopen{}\left\{a\right\}\mathclose{}}{\operatorname{exp}\mathopen{}\left\{b\right\}\mathclose{}}\]

(recall from Algebra 2)

Therefore, according to this model, differences of \(\delta\) in covariate \(x_j\) correspond to rate ratios of \(\operatorname{exp}\mathopen{}\left\{\beta_j \cdot \delta\right\}\mathclose{}\).

That is, letting \(\tilde{X}_{-j}\) denote vector \(\tilde{X}\) with element \(j\) removed:

\[ \begin{aligned} &\phantom{={}} \operatorname{log}\mathopen{}\left\{\operatorname{E}\mathopen{}\left[Y |{\color{red}{X_j = a}}, \tilde{X}_{-j}=\tilde{x}_{-j},T=t\right]\mathclose{}\right\}\mathclose{} \\ &\phantom{={}} -\operatorname{log}\mathopen{}\left\{\operatorname{E}\mathopen{}\left[Y |{\color{red}{X_j = b}}, \tilde{X}_{-j}=\tilde{x}_{-j},T=t\right]\mathclose{}\right\}\mathclose{} \\ &= \operatorname{log}\mathopen{}\left\{t\right\}\mathclose{} + \beta_0 + \beta_1 x_1 + \ldots+ {\color{red}{\beta_j (a)}} + \ldots+ \beta_p x_p \\ &\phantom{={}} -\operatorname{log}\mathopen{}\left\{t\right\}\mathclose{} - \beta_0 - \beta_1 x_1 - \ldots- {\color{red}{\beta_j (b)}} - \ldots- \beta_p x_p \\ &= {\color{red}{\beta_j(a-b)}} \end{aligned} \]

And accordingly,

\[ \begin{aligned} \frac {\operatorname{E}\mathopen{}\left[Y |{\color{red}{X_j = a}}, \tilde{X}_{-j} = \tilde{x}_{-j}, T = t\right]\mathclose{} } { \operatorname{E}\mathopen{}\left[Y |{\color{red}{X_j = b}}, \tilde{X}_{-j}=\tilde{x}_{-j},T=t\right]\mathclose{} } = \operatorname{exp}\mathopen{}\left\{{\color{red}{\beta_j(a-b)}}\right\}\mathclose{} \end{aligned} \]

3 Example: needle-sharing

(adapted from Vittinghoff et al. (2012), §8)

library(tidyverse)
library(haven)
needles <-
  here::here("inst/extdata/needle_sharing.dta") |>
  read_dta() |>
  as_tibble() |>
  mutate(
    hivstat = hivstat |>
      case_match(1 ~ "HIV+", 0 ~ "HIV-") |>
      factor() |>
      relevel(ref = "HIV-"),
    polydrug = polydrug |>
      case_match(1 ~ "multiple drugs used", 0 ~ "one drug used") |>
      factor() |>
      relevel(ref = "one drug used"),
    homeless = homeless |>
      case_match(1 ~ "homeless", 0 ~ "not homeless") |>
      factor() |>
      relevel(ref = "not homeless"),
    sexabuse = sexabuse |>
      case_match(1 ~ "yes", 0 ~ "no") |>
      factor() |>
      relevel(ref = "no"),
    ethn = ethn |>
      factor() |>
      forcats::fct_lump_min(min = 5, other_level = "Other") |>
      relevel(ref = "White"),
    sex = sex |> factor() |> relevel(ref = "M")
  ) |>
  labelled::set_variable_labels(
    "id" = "Participant identifier",
    "sex" = "Biological sex (reference = Male)",
    "ethn" = "Self-reported ethnicity (reference = White; n<5 lumped to Other)",
    "age" = "Age at first interview (years)",
    "dprsn_dx" = "Depression diagnosis indicator",
    "sexabuse" = "History of sexual abuse (1 = yes, 0 = no)",
    "shared_syr" = "Number of shared syringes (in 30 days)",
    "hivstat" = "HIV serostatus (reference = HIV-)",
    "hplsns" = "Hopelessness score",
    "nivdu" = "Number of injections (in 30 days)",
    "shsyryn" = "Shared syringe? (1 = yes, 0 = no)",
    "sqrtnivd" = "sqrt(No. of injections in 30 days)",
    "logshsyr" = "log(No. of shared needles)",
    "polydrug" = "Multiple non-injection drugs used?",
    "sqrtninj" = "sqrt(No. of injections in 30 days); duplicate of sqrtnivd",
    "homeless" = "Homeless? (reference = not homeless)",
    "shsyr" = "Shared needles (30 days); sparse duplicate of shared_syr"
  )

3.1 Introduction

This example, from Vittinghoff et al. (2012, sec. 8.3.1), models risky drug-use behavior — needle sharing — among injection drug users (IDUs), a behavior of public-health interest because shared syringes can transmit blood-borne infections. We examine how needle-sharing frequency relates to participants’ age, sex, ethnicity, housing status (homelessness), drug-use patterns, HIV serostatus, sexual-abuse history, and hopelessness score; injection frequency reflects how often each participant has the opportunity to share a syringe.

These data are a sample of injection drug users distributed with Vittinghoff et al. (2012, sec. 8.3.1) via the UCSF companion website; the original study is not documented in the textbook or the companion materials. One possible (though unconfirmed) source is Tsui et al. (2009), a study of injection-drug-using participants in San Francisco co-authored by one of the textbook’s authors (E. Vittinghoff). There are \(n = 128\) participants and 17 variables. Our outcome is shared_syr, the number of times a participant shared a syringe in the past 30 days.

dict <- tibble(
  variable = names(needles),
  description = labelled::get_variable_labels(needles) |> 
    sapply(function(x) ifelse(is.null(x), "", x)),
)
dict |> pander::pander()

Table 1: Data dictionary for needles data

variable	description
id	Participant identifier
sex	Biological sex (reference = Male)
ethn	Self-reported ethnicity (reference = White; n<5 lumped to Other)
age	Age at first interview (years)
dprsn_dx	Depression diagnosis indicator
sexabuse	History of sexual abuse (1 = yes, 0 = no)
shared_syr	Number of shared syringes (in 30 days)
hivstat	HIV serostatus (reference = HIV-)
hplsns	Hopelessness score
nivdu	Number of injections (in 30 days)
shsyryn	Shared syringe? (1 = yes, 0 = no)
sqrtnivd	sqrt(No. of injections in 30 days)
logshsyr	log(No. of shared needles)
polydrug	Multiple non-injection drugs used?
sqrtninj	sqrt(No. of injections in 30 days); duplicate of sqrtnivd
homeless	Homeless? (reference = not homeless)
shsyr	Shared needles (30 days); sparse duplicate of shared_syr

Table 2 summarizes the demographic and behavioural variables for the participants.

needles |>
  dplyr::select(-id) |>
  gtsummary::tbl_summary()

Table 2: Demographics and behavioural variables for the needle-sharing dataset (all participants).

Characteristic	N = 1281
Biological sex (reference = Male)
M	97 (76%)
F	30 (23%)
Trans	1 (0.8%)
Self-reported ethnicity (reference = White; n<5 lumped to Other)
White	76 (59%)
AA	36 (28%)
Hispanic	10 (7.8%)
Other	6 (4.7%)
Age at first interview (years)	41 (35, 48)
Depression diagnosis indicator
1	88 (70%)
5	38 (30%)
Unknown	2
History of sexual abuse (1 = yes, 0 = no)	17 (14%)
Unknown	4
Number of shared syringes (in 30 days)	0.0 (0.0, 0.0)
Unknown	5
HIV serostatus (reference = HIV-)
HIV-	103 (92%)
HIV+	9 (8.0%)
Unknown	16
Hopelessness score	5.0 (2.0, 10.0)
Number of injections (in 30 days)	90 (40, 120)
Unknown	1
Shared syringe? (1 = yes, 0 = no)	32 (25%)
sqrt(No. of injections in 30 days)	9.5 (6.3, 11.0)
Unknown	1
log(No. of shared needles)	1.61 (0.69, 2.30)
Unknown	101
Multiple non-injection drugs used?
one drug used	109 (85%)
multiple drugs used	19 (15%)
sqrt(No. of injections in 30 days); duplicate of sqrtnivd	9.5 (6.3, 11.0)
Unknown	1
Homeless? (reference = not homeless)
not homeless	63 (51%)
homeless	61 (49%)
Unknown	4
Shared needles (30 days); sparse duplicate of shared_syr	5 (2, 10)
Unknown	101
1 n (%); Median (Q1, Q3)

Before modeling, we sketch a hypothesized causal structure for these variables.

library(dagitty)
library(ggdag)
library(dplyr)

# This is a cross-sectional study: every covariate is measured at a single
# interview, so we cannot causally order the time-varying covariates among
# themselves. We draw their mutual associations -- from shared, unobserved
# earlier circumstances -- as non-directional (bidirected) edges, instead of
# giving each covariate its own separate latent "past state" node.
demographics <- c("age", "sex", "ethn") # fixed at (or determined by) birth
covariates <- c(
  "homeless", "polydrug", "hivstat", "sexabuse", "dprsn_dx", "hplsns"
)

# every unordered pair of time-varying covariates gets a bidirected edge
covariate_pairs <- utils::combn(covariates, 2)
bidirected_edges <- paste(covariate_pairs[1, ], "<->", covariate_pairs[2, ])
demographic_edges <- as.vector(
  outer(demographics, covariates, function(a, b) paste(a, "->", b))
)

needles_dag <- dagitty(paste0("dag {", paste(c(
  # birth-fixed demographics cause the current states, injection frequency,
  # and the outcome
  demographic_edges,
  paste(demographics, "-> nivdu"),
  paste(demographics, "-> shared_syr"),
  # time-varying covariates are mutually associated, with no causal ordering
  bidirected_edges,
  # each time-varying covariate can affect injection frequency and the outcome
  paste(covariates, "-> nivdu"),
  paste(covariates, "-> shared_syr"),
  # injecting is the only mediator: sharing a syringe requires injecting
  "nivdu -> shared_syr",
  # study inclusion is conditioned on; several covariates affect being sampled,
  # so conditioning on `selected` induces associations (selection bias)
  paste(c(demographics, "homeless", "polydrug", "hivstat"), "-> selected")
), collapse = "; "), "}"))

# layout: time flows left -> right (birth-fixed demographics -> current
# covariates -> injection frequency -> shared syringes)
y_cov <- seq(4, -4, length.out = length(covariates))
coordinates(needles_dag) <- list(
  x = c(
    age = -1, sex = -1, ethn = -1,
    setNames(rep(2, length(covariates)), covariates),
    nivdu = 5, shared_syr = 7, selected = 5
  ),
  y = c(
    age = 2, sex = 0, ethn = -2,
    setNames(y_cov, covariates),
    nivdu = 3.2, shared_syr = 0, selected = -5.5
  )
)

needles_tidy <- needles_dag |>
  tidy_dagitty() |>
  mutate(role = case_when(
    name == "homeless" ~ "exposure",
    name == "shared_syr" ~ "outcome",
    name == "nivdu" ~ "mediator",
    name == "selected" ~ "study inclusion (conditioned on)",
    name %in% demographics ~ "demographic (birth-fixed)",
    TRUE ~ "time-varying covariate (current/recent status)"
  ))

# split edges so associations and causal arrows can be styled separately
directed_edges <- function(d) dplyr::filter(d, .data$direction == "->")
association_edges <- function(d) dplyr::filter(d, .data$direction == "<->")

needles_tidy |>
  ggplot() +
  aes(x = x, y = y, xend = xend, yend = yend) +
  # non-directional associations among time-varying covariates: dashed grey arcs
  geom_dag_edges_arc(
    data = association_edges, curvature = 0.4, arrow = NULL,
    edge_colour = "grey65", edge_linetype = "dashed", edge_width = 0.4,
    start_cap = ggraph::circle(5, "mm"), end_cap = ggraph::circle(5, "mm")
  ) +
  # causal (directed) edges: solid arrows
  geom_dag_edges_link(
    data = directed_edges, edge_width = 0.35, edge_colour = "grey20",
    start_cap = ggraph::circle(5, "mm"), end_cap = ggraph::circle(5, "mm")
  ) +
  geom_dag_point(aes(colour = role, shape = role), size = 7) +
  geom_dag_text(colour = "black", size = 1.7) +
  scale_shape_manual(values = c(
    "exposure" = 16, "outcome" = 16, "mediator" = 16,
    "study inclusion (conditioned on)" = 15,
    "demographic (birth-fixed)" = 16,
    "time-varying covariate (current/recent status)" = 16
  )) +
  scale_colour_manual(values = c(
    "exposure" = "#d73027", "outcome" = "#4575b4", "mediator" = "#1a9850",
    "study inclusion (conditioned on)" = "grey30",
    "demographic (birth-fixed)" = "grey55",
    "time-varying covariate (current/recent status)" = "grey80"
  )) +
  guides(colour = guide_legend(nrow = 2), shape = guide_legend(nrow = 2)) +
  theme_dag() +
  theme(legend.position = "bottom", legend.title = element_blank())

Figure 1: Hypothesized causal structure for the needle-sharing example, drawn for a cross-sectional study. The exposure of interest is homelessness (homeless, red) and the outcome is the count of shared syringes (shared_syr, blue). Injection frequency (nivdu, green) is the one genuine mediator: sharing a syringe requires injecting. Demographic variables (age, sex, ethn, darker grey) are fixed at (or determined by) birth, so they sit upstream of every other variable, with arrows into the time-varying covariates, injection frequency, and the outcome. Because every other covariate is measured at the same single interview, we cannot order the time-varying covariates (lighter grey) causally among themselves; their mutual associations — from shared, unobserved earlier circumstances — are drawn as dashed, non-directional edges. Each time-varying covariate may still affect injection frequency and the 30-day outcome, shown as solid arrows into nivdu and shared_syr. The study-inclusion indicator (selected, square) is conditioned on by design; because several covariates influence who is sampled, conditioning on it can induce associations among them (selection bias).

We treat homelessness as the exposure of interest and ask how it affects the number of shared syringes — following Vittinghoff et al. (2012, sec. 8.3.1), who model this outcome on homeless.

Injection frequency nivdu is a mediator: a shared syringe is necessarily one of the injections, so homelessness can raise sharing partly by raising how often a person injects. Because nivdu lies on the causal path from homeless to shared_syr, we do not adjust for it (neither as a covariate nor as an offset); doing so would block that path and estimate only a partial effect. We want the total effect of homelessness on needle-sharing.

Because the study is cross-sectional, the remaining covariates are measured at the same interview as the exposure, so we cannot order them causally relative to homeless or to one another. We draw their mutual associations as non-directional (bidirected) edges and adjust for the measured covariates to control the confounding they carry. A single fitted coefficient is therefore an adjusted association, not an isolated causal effect.

Conditioning on study inclusion (selected) is unavoidable — we only observe people who were sampled — but because several covariates influence who is sampled, it can induce selection bias that no covariate adjustment removes.

needles

Table 3: Needle-sharing data

Figure 2 plots the shared-syringe counts against age, split by the covariate subgroups (point size shows injection frequency).

library(ggplot2)

needles |>
  ggplot() +
  aes(
    x = age,
    y = shared_syr,
    shape = sex,
    col = hivstat
  ) +
  geom_point(aes(size = nivdu), alpha = .5) +
  scale_size_area(max_size = 4) +
  scale_y_continuous(
    transform = "pseudo_log",
    breaks = c(0, 1, 3, 10, 30, 60)
  ) +
  facet_grid(cols = vars(polydrug), rows = vars(homeless)) +
  labs(y = "Number of shared syringes (30 days)") +
  theme(legend.position = "bottom")

Figure 2: Needle-sharing counts by age

Figure 3 shows the underlying counts: the number of shared syringes against the number of injections, with a dashed \(y = x\) line marking the ceiling where every injection involves a shared syringe.

library(ggplot2)

needles |>
  ggplot() +
  aes(
    x = nivdu,
    y = shared_syr,
    shape = sex,
    col = hivstat
  ) +
  geom_point(alpha = .5) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", alpha = .5) +
  # both axes are counts, so pseudo-log both: a shared transform keeps the
  # y = x reference line straight (data diagonal -> panel diagonal)
  scale_x_continuous(
    transform = "pseudo_log",
    breaks = c(0, 3, 10, 30, 100, 300, 900)
  ) +
  scale_y_continuous(
    transform = "pseudo_log",
    breaks = c(0, 1, 3, 10, 30, 60)
  ) +
  facet_grid(cols = vars(polydrug), rows = vars(homeless)) +
  labs(
    x = "Number of injections (in 30 days)",
    y = "Number of shared syringes (in 30 days)"
  ) +
  theme(legend.position = "bottom")

Figure 3: Number of shared syringes versus number of injections (raw counts, 30-day window). The dashed line marks \(y = x\) (every injection involves a shared syringe).

3.2 Covariate counts

Table 4: Counts in needles by sex, housing status, and multiple drug use

needles |>
  dplyr::select(sex, homeless, polydrug) |>
  summary()
#>     sex             homeless                 polydrug  
#>  M    :97   not homeless:63   one drug used      :109  
#>  F    :30   homeless    :61   multiple drugs used: 19  
#>  Trans: 1   NAs         : 4

There’s only one individual with sex = Trans, which unfortunately isn’t enough data to analyze. We will remove that individual:

needles <- needles |> filter(sex != "Trans")

3.3 Model

Following Vittinghoff et al. (2012, sec. 8.3.1), we take homelessness as the exposure of interest and ask how it affects the count of shared syringes, using a Poisson regression (log link) rather than modeling a rate. The remaining covariates enter as adjustments for confounding (Figure 1): because the study is cross-sectional, they are mutually associated with homelessness through shared, unobserved earlier circumstances, which we cannot resolve into a causal order.

Exercise 4 A Poisson rate model would add offset(log(nivdu)), forcing the expected count to be proportional to injection frequency (the coefficient of log(nivdu) fixed at exactly 1). Why not use nivdu as an offset here?

Solution 4. We avoid an offset here for two reasons:

• nivdu is a mediator, not an exposure window. An offset is appropriate when its variable is a known denominator — person-time or population at risk — that is external to the causal question. Injection frequency is instead on the causal path from homelessness to sharing (Figure 1): being homeless can raise sharing partly by raising how often a person injects. Conditioning on it (whether as an offset or as a covariate) would block that path and estimate only a partial effect. We want the total effect of homelessness on needle-sharing.
• It imposes an untested proportionality assumption. Fixing the log(nivdu) coefficient at 1 assumes sharing scales one-for-one with injecting, which the data need not support.

Vittinghoff’s analysis uses neither an offset nor nivdu as a covariate.

Rule of thumb

An offset’s exposure is usually something fixed by the study design — person-time of follow-up, population at risk, area surveyed, number of trials — a known denominator that is exogenous to the question (not caused by the predictors of interest). When the candidate denominator is instead something the participants did, and the predictors also influence it (here, how often they inject), it is a mediator rather than an exposure window, and we model the count directly instead of dividing by it.

We adjust for the measured covariates to control the confounding they carry in Figure 1 — demographic (age, sex, ethn), social (polydrug), serostatus (hivstat), and trauma/mental-health (sexabuse, hopelessness score hplsns) — using main effects only. We exclude the mediator nivdu (above) and the depression indicator dprsn_dx, whose 1/5 coding is undocumented in the source data and so cannot be interpreted.

glm1 <- glm(
  formula = shared_syr ~ homeless +
    age + sex + ethn + polydrug + hivstat + sexabuse + hplsns,
  family = stats::poisson,
  data = needles
)
library(equatiomatic)
equatiomatic::extract_eq(glm1)

\[ \log ({ E( \operatorname{shared_syr} ) }) = \alpha + \beta_{1}(\operatorname{homeless}) + \beta_{2}(\operatorname{age}) + \beta_{3}(\operatorname{sex}_{\operatorname{F}}) + \beta_{4}(\operatorname{ethn}_{\operatorname{AA}}) + \beta_{5}(\operatorname{ethn}_{\operatorname{Hispanic}}) + \beta_{6}(\operatorname{ethn}_{\operatorname{Other}}) + \beta_{7}(\operatorname{polydrug}) + \beta_{8}(\operatorname{hivstat}_{\operatorname{HIV+}}) + \beta_{9}(\operatorname{sex}_{\operatorname{abuseyes}} \times \operatorname{sexabuse}_{\operatorname{yes}}) + \beta_{10}(\operatorname{hplsns}) \]

library(parameters)
glm1 |>
  parameters(exponentiate = TRUE) |>
  print_md()

Table 5: Poisson count model for needle-sharing data

Parameter	IRR	SE	95% CI	z	p
(Intercept)	3.69	1.35	(1.78, 7.46)	3.58	< .001
homeless (homeless)	6.04	1.01	(4.38, 8.45)	10.73	< .001
age	0.96	8.79e-03	(0.94, 0.97)	-4.75	< .001
sex (F)	2.42	0.36	(1.81, 3.23)	5.97	< .001
ethn (AA)	2.59	0.48	(1.80, 3.74)	5.13	< .001
ethn (Hispanic)	1.78	0.57	(0.91, 3.20)	1.81	0.070
ethn (Other)	1.04	0.41	(0.44, 2.08)	0.09	0.927
polydrug (multiple drugs used)	0.19	0.06	(0.10, 0.33)	-5.48	< .001
hivstat (HIV+)	0.11	0.05	(0.04, 0.25)	-4.72	< .001
sexabuse (yes)	0.23	0.07	(0.12, 0.40)	-4.71	< .001
hplsns	0.97	0.01	(0.94, 1.00)	-2.21	0.027

library(ggplot2)

# Predicted expected counts across age for the two homelessness groups,
# holding the other covariates at their reference levels (hplsns at its
# median). The model has no offset, so predict(type = "response") returns
# the expected number of shared syringes directly.
ref_lvl <- function(v) levels(glm1$model[[v]])[1]
pred_grid <- expand.grid(
  age = seq(min(glm1$model$age), max(glm1$model$age), length.out = 100),
  homeless = levels(glm1$model$homeless),
  sex = ref_lvl("sex"),
  ethn = ref_lvl("ethn"),
  polydrug = ref_lvl("polydrug"),
  hivstat = ref_lvl("hivstat"),
  sexabuse = ref_lvl("sexabuse"),
  hplsns = median(glm1$model$hplsns)
)
pred_grid$shared_syr <- predict(glm1, newdata = pred_grid, type = "response")

ggplot(glm1$model) +
  aes(x = age, y = shared_syr, colour = homeless) +
  geom_point(alpha = 0.4) +
  geom_line(data = pred_grid, linewidth = 1) +
  scale_y_continuous(
    transform = "pseudo_log",
    breaks = c(0, 1, 3, 10, 30, 60)
  ) +
  labs(
    y = "Number of shared syringes (30 days)",
    colour = "Housing status"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

Figure 4: Fitted Poisson model: expected sharing by age and housing status

3.4 Interpreting the coefficients

For a Poisson model with a log link, exponentiating a coefficient turns it into a ratio of expected counts (an incidence rate ratio, IRR, in Vittinghoff et al. (2012)’s terminology): \(\operatorname{exp}\mathopen{}\left\{\beta_j\right\}\mathclose{}\) is the multiplicative change in the expected number of shared syringes per one-unit increase in \(x_j\), holding the other covariates fixed. Our focal estimate is the IRR for homelessness; the remaining coefficients are adjustments for confounding, not effects of separate interest. Because we do not condition on the mediator nivdu, the homelessness IRR is a total effect. The fitted ratios are reported in Table 5.

irr <- exp(coef(glm1))

Baseline (intercept). \(\operatorname{exp}\mathopen{}\left\{\hat\beta_0\right\}\mathclose{} = 3.69\) is the expected count of shared syringes for the reference participant (male, not homeless, White, one drug used, HIV-negative, no sexual-abuse history, hopelessness score 0), extrapolated to age 0. Age 0 is outside the data, so the intercept is a mathematical anchor rather than a directly interpretable quantity.

Age. \(\operatorname{exp}\mathopen{}\left\{\hat\beta_{\text{age}}\right\}\mathclose{} = 0.957\) per year: essentially no association between age and the expected count.

Sex. \(\operatorname{exp}\mathopen{}\left\{\hat\beta_{\text{sexF}}\right\}\mathclose{} = 2.42\): women’s expected count of shared syringes is about 2.42 times men’s, adjusting for the other covariates.

Homelessness (the focal exposure). \(\operatorname{exp}\mathopen{}\left\{\hat\beta_{\text{homeless}}\right\}\mathclose{} = 6.04\): homeless participants share syringes about 6.04 times as often as the non-homeless, adjusting for the other covariates. Vittinghoff et al. (2012, sec. 8.3.1) reports an unadjusted IRR of about 3.3 for homelessness in a single-predictor model; adjusting for the confounder-proxies here moves the estimate to about 6.

Multiple-drug use. \(\operatorname{exp}\mathopen{}\left\{\hat\beta_{\text{polydrug}}\right\}\mathclose{} = 0.19\): multiple-drug users have a lower expected count than single-drug users, though only 19 participants used multiple drugs, so this estimate is imprecise.

HIV serostatus. \(\operatorname{exp}\mathopen{}\left\{\hat\beta_{\text{hivstatHIV+}}\right\}\mathclose{} = 0.11\): HIV-positive participants have a lower expected count of shared syringes than HIV-negative participants, adjusting for the other covariates. Only 8 participants are HIV+ in this sample, so this estimate is imprecise.

Ethnicity. The ethn coefficients compare each retained group to the White reference (small groups were lumped into “Other”); they adjust for any confounding of the homelessness effect by ethnicity but are not of separate interest.

Sexual-abuse history. \(\operatorname{exp}\mathopen{}\left\{\hat\beta_{\text{sexabuse}}\right\}\mathclose{} = 0.23\): participants with a history of sexual abuse share syringes about 0.23 times as often as those without, adjusting for the other covariates.

Hopelessness. \(\operatorname{exp}\mathopen{}\left\{\hat\beta_{\text{hplsns}}\right\}\mathclose{} = 0.97\) per point: each one-point increase in the hopelessness score multiplies the expected count by 0.97, adjusting for the other covariates.

Overdispersion caveat

These data are highly overdispersed — the variance far exceeds the mean — so the model-based Poisson standard errors are too small. Vittinghoff et al. (2012, sec. 8.3.1) recommends robust standard errors or a negative-binomial model for this example.

4 Inference for count regression models

Confidence intervals for regression coefficients and rate ratios

A Wald 95% confidence interval for a single coefficient \(\beta_j\) is:

\[ \beta_j \in \mathopen{}\left[\hat\beta_j \pm z_{1 - \frac{\alpha}{2}}\cdot \mathop{\widehat{\operatorname{se}}}\nolimits\mathopen{}\left(\hat\beta_j\right)\mathclose{}\right]\mathclose{} \]

where \(z_{1-\alpha/2} \approx 1.96\) for \(\alpha = 0.05\).

Because the log-rate scale is related to the rate scale by exponentiation, we obtain a confidence interval for the rate ratio \(e^{\beta_j}\) by exponentiating both endpoints:

\[ e^{\beta_j} \in \mathopen{}\left[ \operatorname{exp}\mathopen{}\left\{\hat\beta_j - z_{1 - \frac{\alpha}{2}}\cdot \mathop{\widehat{\operatorname{se}}}\nolimits\mathopen{}\left(\hat\beta_j\right)\mathclose{}\right\}\mathclose{},\; \operatorname{exp}\mathopen{}\left\{\hat\beta_j + z_{1 - \frac{\alpha}{2}}\cdot \mathop{\widehat{\operatorname{se}}}\nolimits\mathopen{}\left(\hat\beta_j\right)\mathclose{}\right\}\mathclose{} \right]\mathclose{} \]

Hypothesis tests for regression coefficients

To test \(H_0: \beta_j = \beta_0\) against a one- or two-sided alternative, compute the Wald \(z\)-statistic:

\[ z = \frac{\hat \beta_j - \beta_0}{\mathop{\widehat{\operatorname{se}}}\nolimits\mathopen{}\left(\hat\beta_j\right)\mathclose{}} \]

and compare \(z\) (one-sided) or \(|z|\) (two-sided) to the tails of the standard Gaussian distribution. The most common null hypothesis is \(H_0: \beta_j = 0\), i.e., that covariate \(j\) has no association with the outcome rate.

Comparing nested models

To compare a smaller model \(M_0\) (with \(p_0\) parameters) to a larger model \(M_1\) (with \(p_1 > p_0\) parameters), use the likelihood ratio test statistic:

\[ G^2 = 2\bigl[\hat\ell_1 - \hat\ell_0\bigr] \]

where \(\hat\ell_1\) and \(\hat\ell_0\) are the maximized log-likelihoods of \(M_1\) and \(M_0\) respectively. (Here the subscripts index the two models; they are unrelated to the scalar null value \(\beta_0\) used in the Wald test above.)

Under \(H_0\) that the additional \(p_1 - p_0\) parameters are all zero, \(G^2 \ \dot{\sim} \ \chi^2_{p_1 - p_0}\).

See Dobson and Barnett (2018, chap. 9) and Vittinghoff et al. (2012, sec. 8.1) for details.

5 Prediction

\[ \begin{aligned} \hat y &\stackrel{\text{def}}{=}\hat{\operatorname{E}}[Y|\tilde{X}= \tilde{x},T=t]\\ &=\hat\mu(\tilde{x}, t)\\ &=\hat\lambda(\tilde{x}) \cdot t\\ &=\operatorname{exp}\mathopen{}\left\{\hat\eta(\tilde{x})\right\}\mathclose{} \cdot t\\ &=\operatorname{exp}\mathopen{}\left\{\tilde{x}'\hat{\tilde{\beta}}\right\}\mathclose{} \cdot t \end{aligned} \]

6 Diagnostics

Residuals

Observation residuals

\[e \stackrel{\text{def}}{=}y - \hat y\]

Pearson residuals

\[r = \frac{e}{\mathop{\widehat{\operatorname{se}}}\nolimits\mathopen{}\left(e\right)\mathclose{}} \approx \frac{e}{\sqrt{\hat y}}\]

Standardized Pearson residuals

\[r_p = \frac{r}{\sqrt{1-h}}\] where \(h\) is the “leverage” (which we will continue to leave undefined).

Deviance residuals

\[ d_k = \operatorname{sign}(y - \hat y)\left\{\sqrt{2[\ell_{\text{full}}(y) - \ell(\hat\beta; y)]}\right\} \]

Note

\[\operatorname{sign}(x) \stackrel{\text{def}}{=}\frac{x}{|x|}\] In other words:

• \(\operatorname{sign}(x) = -1\) if \(x < 0\)
• \(\operatorname{sign}(x) = 0\) if \(x = 0\)
• \(\operatorname{sign}(x) = 1\) if \(x > 0\)

library(ggfortify)
autoplot(glm1)

Table 6: Diagnostics for Poisson model

7 Zero-inflation

Models for zero-inflated counts

We assume a latent (unobserved) binary variable, \(Z\), which we model using logistic regression:

\[P(Z=1|X=x) = \pi(x) = \operatorname{expit}(\gamma_0 + \gamma_1 x_1 +...)\]

According to this model, if \(Z=1\), then \(Y\) will always be zero, regardless of \(X\) and \(T\):

\[P(Y=0|Z=1,X=x,T=t) = 1\]

Otherwise (if \(Z=0\)), \(Y\) will have a Poisson distribution, conditional on \(X\) and \(T\), as above.

Even though we never observe \(Z\), we can estimate the parameters \(\gamma_0\)-\(\gamma_p\), via maximum likelihood:

\[ \begin{aligned} \operatorname{P}(Y=y|X=x,T=t) &= \operatorname{P}(Y=y,Z=1|...) + \operatorname{P}(Y=y,Z=0|...) \end{aligned} \] (by the Law of Total Probability)

where \[ \begin{aligned} P(Y=y,Z=z|...) &= P(Y=y|Z=z,...)P(Z=z|...) \end{aligned} \]

Exercise 5 Expand \(P(Y=0|X=x,T=t)\), \(P(Y=1|X=x,T=t)\) and \(P(Y=y|X=x,T=t)\) into expressions involving \(P(Z=1|X=x)\) and \(P(Y=y|Z=0,X=x,T=t)\).

Solution. Let \(\pi = \operatorname{P}(Z=1|X=x)\) and \(\mu_0 = \operatorname{E}[Y|Z=0,X=x,T=t]\).

\(P(Y=0)\): \(Y=0\) occurs either because \(Z=1\) (always zero) or because \(Z=0\) and the Poisson draw equals 0:

\[ \begin{aligned} \operatorname{P}(Y=0|X=x,T=t) &= \operatorname{P}(Z=1|X=x) + \operatorname{P}(Z=0|X=x)\,\operatorname{P}(Y=0|Z=0,X=x,T=t)\\ &= \pi + (1-\pi)\,e^{-\mu_0} \end{aligned} \]

\(P(Y=1)\): \(Z=1\) can never produce \(Y=1\), so:

\[ \operatorname{P}(Y=1|X=x,T=t) = (1-\pi)\,\operatorname{P}(Y=1|Z=0,X=x,T=t) = (1-\pi)\,\mu_0 e^{-\mu_0} \]

\(P(Y=y)\) for \(y \geq 1\): Identical reasoning gives

\[ \operatorname{P}(Y=y|X=x,T=t) = (1-\pi)\,\frac{\mu_0^y e^{-\mu_0}}{y!} \]

Exercise 6 Derive the expected value and variance of \(Y\), conditional on \(X\) and \(T\), as functions of \(\pi = P(Z=1|X=x)\) and \(\mu_0 = \operatorname{E}[Y|Z=0,X=x,T=t]\).

Solution. Let \(\pi = \operatorname{P}(Z=1|X=x)\) and \(\mu_0 = \operatorname{E}[Y|Z=0,X=x,T=t]\).

Expected value. By the Law of Total Expectation (conditioning on \(Z\), within the subpopulation \(\{X=x, T=t\}\)):

\[ \begin{aligned} \operatorname{E}[Y|X=x,T=t] &= \operatorname{E}[Y|Z=1,X=x,T=t]\,\pi + \operatorname{E}[Y|Z=0,X=x,T=t]\,(1-\pi)\\ &= 0 \cdot \pi + \mu_0(1-\pi)\\ &= (1-\pi)\,\mu_0 \end{aligned} \]

The substitution \(\operatorname{E}[Y|Z=0,X=x,T=t] = \mu_0\) follows immediately from the definition of \(\mu_0\) above.

Variance. By the Law of Total Variance. To reduce clutter, we suppress the \((X=x, T=t)\) conditioning in the intermediate steps below: every expectation and variance is taken within the subpopulation \(\{X=x, T=t\}\), and we restore the explicit conditioning in the final line.

\[ \operatorname{Var}\mathopen{}\left(Y\right)\mathclose{} = \operatorname{E}[\operatorname{Var}\mathopen{}\left(Y|Z\right)\mathclose{}] + \operatorname{Var}\mathopen{}\left(\operatorname{E}[Y|Z]\right)\mathclose{} \]

For the first term, since \(\operatorname{Var}\mathopen{}\left(Y|Z=1\right)\mathclose{}=0\) and \(\operatorname{Var}\mathopen{}\left(Y|Z=0\right)\mathclose{}=\mu_0\) (Poisson):

\[ \operatorname{E}[\operatorname{Var}\mathopen{}\left(Y|Z\right)\mathclose{}] = 0 \cdot \pi + \mu_0(1-\pi) = (1-\pi)\mu_0 \]

For the second term, \(\operatorname{E}[Y|Z]\) takes the value 0 (with prob \(\pi\)) or \(\mu_0\) (with prob \(1-\pi\)), so:

\[ \begin{aligned} \operatorname{Var}\mathopen{}\left(\operatorname{E}[Y|Z]\right)\mathclose{} &= \pi(0 - (1-\pi)\mu_0)^2 + (1-\pi)(\mu_0 - (1-\pi)\mu_0)^2\\ &= \pi(1-\pi)^2\mu_0^2 + (1-\pi)\pi^2\mu_0^2\\ &= \pi(1-\pi)\mu_0^2[\,(1-\pi)+\pi\,]\\ &= \pi(1-\pi)\mu_0^2 \end{aligned} \]

Combining:

\[ \operatorname{Var}\mathopen{}\left(Y|X=x,T=t\right)\mathclose{} = (1-\pi)\mu_0 + \pi(1-\pi)\mu_0^2 = (1-\pi)\mu_0\bigl(1 + \pi\mu_0\bigr) \]

Since \((1-\pi)\mu_0\bigl(1+\pi\mu_0\bigr) \geq (1-\pi)\mu_0 = \operatorname{E}[Y|X=x,T=t]\) for any \(\pi > 0\), zero-inflated count models always exhibit overdispersion relative to a Poisson model with the same mean.

8 Over-dispersion

The Poisson distribution model forces the variance to equal the mean. In practice, many count distributions will have a variance substantially larger than the mean (or occasionally, smaller).

Definition 1 (Overdispersion) A random variable \(X\) is overdispersed relative to a model \(\operatorname{p}(X=x)\) if its empirical variance in a dataset is larger than the value predicted by the fitted model \(\hat{\operatorname{p}}(X=x)\).

c.f. Dobson and Barnett (2018) §3.2.1, 7.7, 9.8; Vittinghoff et al. (2012) §8.1.5; and https://en.wikipedia.org/wiki/Overdispersion.

When we encounter overdispersion, we can try to reduce the residual variance by adding more covariates.

Note

Logistic regression is named after the (inverse) link function. Poisson regression is named after the outcome distribution. I think this naming convention reflects the strongest (most questionable assumption) in the model. In binary data regression, the outcome distribution essentially must be Bernoulli (or Binomial), but the link function could be logit, log, identity, probit, or something more unusual. In count data regression, the outcome distribution could have many different shapes, but the link function will probably end up being log, so that covariates have multiplicative effects on the rate.

8.1 Negative binomial models

There are alternatives to the Poisson model. Most notably, the negative binomial model.

The negative binomial distribution is a common alternative to the Poisson distribution for count outcomes. It adds a dispersion parameter that allows the variance to exceed the mean, making it more flexible when overdispersion is present. We can still model \(\mu\) as a function of \(X\) and \(T\) as before, and we can combine this model with zero-inflation (as the conditional distribution for the non-zero component).

Example: needle-sharing

library(MASS) # for glm.nb()
glm1_nb <- glm.nb(
  formula = shared_syr ~ homeless +
    age + sex + ethn + polydrug + hivstat + sexabuse + hplsns,
  data = needles
)

equatiomatic::extract_eq(glm1_nb)

\[ \log ({ E( \operatorname{shared_syr} ) }) = \alpha + \beta_{1}(\operatorname{homeless}) + \beta_{2}(\operatorname{age}) + \beta_{3}(\operatorname{sex}_{\operatorname{F}}) + \beta_{4}(\operatorname{ethn}_{\operatorname{AA}}) + \beta_{5}(\operatorname{ethn}_{\operatorname{Hispanic}}) + \beta_{6}(\operatorname{ethn}_{\operatorname{Other}}) + \beta_{7}(\operatorname{polydrug}) + \beta_{8}(\operatorname{hivstat}_{\operatorname{HIV+}}) + \beta_{9}(\operatorname{sex}_{\operatorname{abuseyes}} \times \operatorname{sexabuse}_{\operatorname{yes}}) + \beta_{10}(\operatorname{hplsns}) \]

Table 7: Negative binomial model for needle-sharing data

summary(glm1_nb)
#> 
#> Call:
#> glm.nb(formula = shared_syr ~ homeless + age + sex + ethn + polydrug + 
#>     hivstat + sexabuse + hplsns, data = needles, init.theta = 0.0903642555, 
#>     link = log)
#> 
#> Coefficients:
#>                             Estimate Std. Error z value Pr(>|z|)   
#> (Intercept)                   2.8948     1.9293    1.50   0.1335   
#> homelesshomeless              2.5534     0.8024    3.18   0.0015 **
#> age                          -0.0729     0.0447   -1.63   0.1027   
#> sexF                         -0.5545     0.8869   -0.63   0.5318   
#> ethnAA                        0.9100     0.9416    0.97   0.3339   
#> ethnHispanic                  2.6222     1.5835    1.66   0.0977 . 
#> ethnOther                    -1.5963     2.1745   -0.73   0.4629   
#> polydrugmultiple drugs used  -3.7062     1.3330   -2.78   0.0054 **
#> hivstatHIV+                  -3.4809     1.5036   -2.32   0.0206 * 
#> sexabuseyes                  -1.5612     1.2484   -1.25   0.2111   
#> hplsns                       -0.0449     0.0685   -0.66   0.5119   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for Negative Binomial(0.0899) family taken to be 1)
#> 
#>     Null deviance: 61.511  on 108  degrees of freedom
#> Residual deviance: 52.284  on  98  degrees of freedom
#>   (18 observations deleted due to missingness)
#> AIC: 261.5
#> 
#> Number of Fisher Scoring iterations: 25
#> 
#> 
#>               Theta:  0.0904 
#>           Std. Err.:  0.0239 
#> Warning while fitting theta: alternation limit reached 
#> 
#>  2 x log-likelihood:  -237.5430

tibble(name = names(coef(glm1)), poisson = coef(glm1), nb = coef(glm1_nb))

Table 8: Poisson versus Negative Binomial Regression coefficient estimates

zero-inflation

library(glmmTMB)
zinf_pois <- glmmTMB(
  family = "poisson",
  data = needles,
  formula = shared_syr ~ homeless +
    age + sex + ethn + polydrug + hivstat + sexabuse + hplsns,
  # the zero-inflation submodel uses the exposure, matching Vittinghoff's
  # `inflate(i.homeless)` (@vittinghoff2e, §8.3.1)
  ziformula = ~ homeless
)

zinf_pois |>
  parameters(exponentiate = TRUE) |>
  print_md()

Table 9: Zero-inflated Poisson model

Fixed Effects

Parameter	IRR	SE	95% CI	z	p
(Intercept)	0.42	0.42	(0.06, 2.94)	-0.87	0.385
homeless (homeless)	4.13	0.90	(2.70, 6.32)	6.54	< .001
age	1.02	0.02	(0.99, 1.06)	1.14	0.256
sex (F)	10.01	3.44	(5.10, 19.63)	6.70	< .001
ethn (AA)	11.21	4.24	(5.35, 23.51)	6.40	< .001
ethn (Hispanic)	0.62	0.25	(0.28, 1.36)	-1.20	0.230
ethn (Other)	1.25	0.64	(0.46, 3.42)	0.43	0.664
polydrug (multiple drugs used)	0.98	0.50	(0.36, 2.68)	-0.03	0.974
hivstat (HIV+)	0.36	0.18	(0.14, 0.97)	-2.02	0.044
sexabuse (yes)	0.09	0.03	(0.05, 0.18)	-6.87	< .001
hplsns	1.05	0.03	(0.99, 1.12)	1.66	0.096

Zero-Inflation

Parameter	Odds Ratio	SE	95% CI	z	p
(Intercept)	4.84	1.92	(2.22, 10.55)	3.97	< .001
homeless (homeless)	0.53	0.27	(0.20, 1.44)	-1.24	0.216

Another R package for zero-inflated models is pscl (Zeileis et al. (2008)).

zero-inflated negative binomial model

zinf_nb <- glmmTMB(
  family = nbinom2,
  data = needles,
  formula = shared_syr ~ homeless +
    age + sex + ethn + polydrug + hivstat + sexabuse + hplsns,
  ziformula = ~ homeless
)

zinf_nb |>
  parameters(exponentiate = TRUE) |>
  print_md()

Table 10: Zero-inflated negative binomial model

Fixed Effects

Parameter	IRR	SE	95% CI	z	p
(Intercept)	0.58	2.93	(2.76e-05, 12031.38)	-0.11	0.914
homeless (homeless)	2.66	1.94	(0.64, 11.08)	1.34	0.180
age	1.02	0.10	(0.85, 1.23)	0.26	0.798
sex (F)	8.11	13.94	(0.28, 235.47)	1.22	0.223
ethn (AA)	9.40	12.29	(0.72, 121.95)	1.71	0.087
ethn (Hispanic)	0.90	1.38	(0.04, 18.04)	-0.07	0.945
ethn (Other)	1.46	2.97	(0.03, 77.95)	0.19	0.851
polydrug (multiple drugs used)	0.79	2.02	(5.15e-03, 120.47)	-0.09	0.926
hivstat (HIV+)	0.42	0.52	(0.04, 4.72)	-0.70	0.481
sexabuse (yes)	0.09	0.09	(0.01, 0.65)	-2.40	0.016
hplsns	1.03	0.12	(0.83, 1.29)	0.30	0.767

Zero-Inflation

Parameter	Odds Ratio	SE	95% CI	z	p
(Intercept)	4.70	1.94	(2.09, 10.56)	3.75	< .001
homeless (homeless)	0.49	0.26	(0.17, 1.40)	-1.33	0.182

Dispersion

Parameter	Coefficient	95% CI
(Intercept)	1.49	(0.54, 4.09)

Both zero-inflated models keep the zero-inflation (structural-zero) submodel parsimonious — just the exposure, ziformula = ~ homeless — matching Vittinghoff et al. (2012)’s inflate(i.homeless) [§8.3.1].

A richer zero-inflation submodel is poorly identified here. With the negative binomial in particular, the overdispersion parameter and the zero-inflation component both compete to explain the excess zeros, so adding several covariates to the submodel makes the logistic fit separate: its coefficients diverge and the exponentiated estimates (odds ratios) run off to \(\pm\infty\). Keeping the submodel small avoids that.

8.2 Quasipoisson

Another way to handle overdispersion — rather than switching to the negative binomial distributional family — is the “quasipoisson” approach. It is a method-of-moments-type inference method: rather than specifying a full probability distribution and fitting by maximum likelihood, it specifies only the mean-variance relationship \(\operatorname{Var}\mathopen{}\left(Y\right)\mathclose{} = \mu\theta\), and estimates \(\theta\) accordingly. This approach is simpler to implement, but provides less information than the full negative binomial likelihood.

See ?quasipoisson in R for more.

9 More on count regression

• https://bookdown.org/roback/bookdown-BeyondMLR/ch-poissonreg.html

Exercises

Exercise 7 (Interpreting Poisson regression coefficients) (adapted from Dunn and Smyth (2018), Chapter 5, and Dobson and Barnett (2018), Chapter 9)

Consider a Poisson log-linear model for the expected number of events \(\mu_i\):

\[ \operatorname{log}\mathopen{}\left\{\mu_i\right\}\mathclose{} = \beta_0 + \beta_1 x_i, \]

where \(x_i\) is a binary indicator (\(x_i = 0\) or \(x_i = 1\)).

(a) Express \(\mu_i\) as a function of \(x_i\).

(b) Interpret \(e^{\beta_0}\).

(c) Interpret \(e^{\beta_1}\).

(d) If \(\hat\beta_0 = 1.2\) and \(\hat\beta_1 = 0.5\), compute the estimated mean event count for \(x_i = 0\) and \(x_i = 1\).

Solution. (a)

\[ \mu_i = e^{\beta_0 + \beta_1 x_i} = e^{\beta_0} \cdot (e^{\beta_1})^{x_i} \]

For \(x_i = 0\): \(\mu_0 = e^{\beta_0}\). For \(x_i = 1\): \(\mu_1 = e^{\beta_0 + \beta_1}\).

(b)

\(e^{\beta_0}\) is the expected mean count when \(x_i = 0\) (the reference group).

(c)

\[ e^{\beta_1} = \frac{\mu_1}{\mu_0} = \frac{e^{\beta_0+\beta_1}}{e^{\beta_0}} \]

\(e^{\beta_1}\) is the rate ratio (or count ratio): the multiplicative factor by which the expected count changes when \(x_i\) increases from 0 to 1.

If \(\beta_1 > 0\), the group with \(x_i = 1\) has a higher expected count.

(d)

For \(x_i = 0\): \[ \hat\mu_0 = e^{1.2} \approx 3.32 \]

For \(x_i = 1\): \[ \hat\mu_1 = e^{1.2 + 0.5} = e^{1.7} \approx 5.47 \]

The estimated rate ratio is \(e^{0.5} \approx 1.65\), meaning the group with \(x_i = 1\) has about 65% more events on average.

Exercise 8 (Score equations for a Poisson GLM) (adapted from Dobson and Barnett (2018), Chapter 4, and Dunn and Smyth (2018), Chapter 4)

Consider the Poisson log-linear model \(\operatorname{log}\mathopen{}\left\{\mu_i\right\}\mathclose{} = \beta_0 + \beta_1 x_i\), with \(Y_i \mid x_i \ \sim_{\perp\!\!\!\perp}\ \operatorname{Pois}(\mu_i)\).

(a) Write the log-likelihood \(\ell(\beta_0, \beta_1; \tilde{y}, \tilde{x})\).

(b) Derive the score equations \(\frac{\partial}{\partial \beta_0}\ell = 0\) and \(\frac{\partial}{\partial \beta_1}\ell = 0\).

(c) Interpret the score equations: what condition on the fitted values \(\hat\mu_i\) do they imply?

Solution. (a)

\[ \begin{aligned} \ell(\beta_0, \beta_1; \tilde{y}, \tilde{x}) &= \sum_{i=1}^n \mathopen{}\left(y_i \operatorname{log}\mathopen{}\left\{\mu_i\right\}\mathclose{} - \mu_i - \operatorname{log}\mathopen{}\left\{y_i!\right\}\mathclose{}\right)\mathclose{} \\ &\propto \sum_{i=1}^n \mathopen{}\left(y_i \operatorname{log}\mathopen{}\left\{\mu_i\right\}\mathclose{} - \mu_i\right)\mathclose{} \end{aligned} \]

where \(\operatorname{log}\mathopen{}\left\{\mu_i\right\}\mathclose{} = \beta_0 + \beta_1 x_i\), so \(\mu_i = e^{\beta_0 + \beta_1 x_i}\).

(b)

Using the chain rule with \(\frac{\partial}{\partial \beta_j}\mu_i = \mu_i x_{ij}\) (where \(x_{i0} = 1\) and \(x_{i1} = x_i\)):

\[ \frac{\partial}{\partial \beta_j}\ell = \sum_{i=1}^n \mathopen{}\left(\frac{y_i}{\mu_i} - 1\right)\mathclose{}\mu_i x_{ij} = \sum_{i=1}^n (y_i - \mu_i) x_{ij} \]

Setting the score equations to zero:

\[ \begin{aligned} \frac{\partial}{\partial \beta_0}\ell = 0 &\;\Rightarrow\; \sum_{i=1}^n (y_i - \mu_i) = 0 \\ \frac{\partial}{\partial \beta_1}\ell = 0 &\;\Rightarrow\; \sum_{i=1}^n x_i(y_i - \mu_i) = 0 \end{aligned} \]

(c)

The first equation says \(\sum_i y_i = \sum_i \hat\mu_i\): the total fitted count equals the total observed count.

The second equation says \(\sum_i x_i y_i = \sum_i x_i \hat\mu_i\): the fitted counts are balanced against observed counts, weighted by \(x_i\).

More generally, these score equations say that the residuals \((y_i - \hat\mu_i)\) are orthogonal to each predictor column: for each predictor \(j\), the residual vector \((\tilde{y}- \hat{\tilde{\mu}})\) satisfies \({\tilde{x}_{(j)}}^{\top}(\tilde{y}- \hat{\tilde{\mu}}) = 0\), where \(\tilde{x}_{(j)} = (x_{1j}, \ldots, x_{nj})\) is the column of \(j\)-th predictor values across observations. This is the GLM analogue of the OLS normal equations.

References

Dobson, Annette J, and Adrian G Barnett. 2018. An Introduction to Generalized Linear Models. 4th ed. CRC press. https://doi.org/10.1201/9781315182780.

Dunn, Peter K, and Gordon K Smyth. 2018. Generalized Linear Models with Examples in R. Vol. 53. Springer. https://link.springer.com/book/10.1007/978-1-4419-0118-7.

Tsui, Judith I., Eric Vittinghoff, Judith A. Hahn, Jennifer L. Evans, Peter J. Davidson, and Kimberly Page. 2009. “Risk Behaviors After Hepatitis C Virus Seroconversion in Young Injection Drug Users in San Francisco.” Drug and Alcohol Dependence 105 (1-2): 160–63. https://doi.org/10.1016/j.drugalcdep.2009.05.022.

Vittinghoff, Eric, David V Glidden, Stephen C Shiboski, and Charles E McCulloch. 2012. Regression Methods in Biostatistics: Linear, Logistic, Survival, and Repeated Measures Models. 2nd ed. Springer. https://doi.org/10.1007/978-1-4614-1353-0.

Zeileis, Achim, Christian Kleiber, and Simon Jackman. 2008. “Regression Models for Count Data in R.” Journal of Statistical Software 27 (8). https://www.jstatsoft.org/v27/i08/.