1 Why Count Models Matter in Political Science

Key Takeaways for Political Scientists

Never use OLS for count data - It violates assumptions and can predict negative counts

Start with Poisson, test for overdispersion - Don’t jump to complex models

Theory determines model choice - ZIP vs Hurdle depends on your causal story

Always use offsets for exposure - Population, time, or space matters

IRRs are your friend - $\exp(\beta)$ gives multiplicative change in expected count

Zero-inflation ≠ Excess zeros - Check specification before adding complexity

Political science runs on counts. International relations scholars count militarized disputes, wars, and treaties. Comparative politics researchers count coups, protests, and government turnovers. Americanists count bill sponsorships, roll call votes, and Supreme Court reversals. These aren’t continuous variables that happen to be integers, they’re fundamentally discrete processes that generate non-negative whole numbers.

1.1 The Problem with Standard Approaches

You’ve been trained in OLS. When you see a numeric dependent variable, your instinct is to run a linear regression. This instinct will fail you with counts, sometimes catastrophically.

1.1.1 Why OLS Fails for Counts

Consider modeling the number of protests in a country-year. The data might look like: 0, 0, 1, 0, 0, 0, 3, 0, 15, 0, 0, 2, 0, 0, 45, 0, 1, 0, 0, 0. Most observations cluster at zero, a few countries have modest counts, and occasionally you get an Arab Spring with dozens of events. OLS assumes your errors are normally distributed around a linear conditional mean. But count data violates every assumption:

Non-negativity constraint: OLS can predict negative protests. This isn’t just silly, it’s impossible.
Discrete outcomes: You can’t have 2.7 coups. Counts are integers.
Heteroskedasticity: Variance increases with the mean (countries with more protests have more variable protest counts).
Skewness: Most observations pile up at low values with a long right tail.

1.1.1.1 The Bias Question: Does OLS Actually Produce Biased Coefficients?

This is where things get subtle and many textbooks oversimplify. The debate isn’t just “OLS bad, Poisson good”, it’s about understanding exactly what goes wrong and when.

The seminal paper by Santos Silva and Tenreyro (2006, Review of Economics and Statistics) exposed a critical flaw: when you log-transform count data and run OLS, you get biased estimates. But why? The problem is that the logarithm is a nonlinear transformation, and Jensen’s inequality tells us that E[log(Y)] ≠ log(E[Y]). When you run OLS on log(Y), you’re modeling E[log(Y)], but what you usually care about is log(E[Y]), the log of the expected count, not the expected log-count. This distinction might seem pedantic, but it has massive consequences for your estimates.

However, recent work has refined this understanding, revealing when OLS might work and when it definitely won’t:

When OLS is (surprisingly) unbiased: Wooldridge (2023, Econometric Analysis of Cross Section and Panel Data, 3rd ed.) demonstrates a counterintuitive result: under certain conditions, OLS coefficients can be unbiased for the average partial effects even with count data. The key word is “average”, OLS gives you the effect on the conditional mean function averaged across the distribution of covariates. But here’s the catch: while the coefficients might be unbiased in this specific sense, the standard errors will be wrong (heteroskedasticity), you lose efficiency (larger standard errors than necessary), and you still can’t interpret the coefficients as you would in a properly specified count model. It’s like using a hammer to drive in a screw, it might work, but why would you?
The functional form problem: Mullahy (2024, Journal of Econometrics) clarifies what’s really going on. The fundamental issue isn’t that Y is a count, it’s that the true conditional mean function is exponential: $E[Y|X] = \exp(X\beta)$. OLS assumes $E[Y|X] = X\beta$, forcing a linear relationship on an exponential curve. This creates specification bias that gets worse as the relationship becomes more nonlinear. Think about it: when protest rates are low (most autocracies), the exponential is nearly linear, so OLS might approximate reality. But when rates are high (during the Arab Spring), the exponential curves sharply upward while OLS plods along linearly, creating massive prediction errors. The bias isn’t constant, it varies across the covariate space, being worst exactly where the most interesting political variation occurs.
Heteroskedasticity-induced bias: Chen, Ponomareva, and Tamer (2023, Quantitative Economics) prove an even more damning result. Count data necessarily exhibits heteroskedasticity, the variance increases with the mean. This isn’t an annoyance we can fix with robust standard errors; it fundamentally breaks OLS. When the variance function depends on covariates (and with counts, $\text{Var}(Y|X) = \mu + \alpha\mu^2$ where $\mu = E[Y|X]$), OLS produces biased estimates of even the linear projection coefficients. Yes, you can use heteroskedasticity-robust standard errors to get valid inference for these biased coefficients, but what’s the point? You’re getting precise estimates of the wrong quantities. It’s like having a very accurate map of the wrong city, technically correct but practically useless for navigation.

1.1.1.2 Connecting to Causal Inference: What Are We Actually Estimating?

Modern causal inference has clarified what we mean by “the effect” of a treatment or covariate. This clarification is crucial for understanding the bias debate with count models.

The causal estimand: In causal inference terms, we typically care about the Average Treatment Effect (ATE) or Average Marginal Effect (AME):

\[ \text{AME} = E\left[\frac{\partial E[Y|X,D]}{\partial D}\right] \]

This is the expected change in the outcome when we change the treatment/covariate, averaged over the population. With count data and an exponential mean function, this becomes:

\[ \text{AME} = E[\beta_D \exp(X\beta)] = \beta_D \times E[\exp(X\beta)] \]

Notice that the marginal effect depends on both $\beta_D$ (the coefficient) AND the expected value of the exponential function. This creates several problems for OLS:

Problem 1: OLS estimates the wrong functional form

OLS with count data estimates something like: \[ \hat{\beta}_{OLS} \approx \arg\min E[(Y - X\beta)^2] \]

This gives you the linear projection of Y onto X. But when the true relationship is $E[Y|X] = \exp(X\beta)$, the linear projection doesn’t equal the marginal effect. The OLS coefficient $\hat{\beta}_{OLS}$ is trying to approximate a curved relationship with a straight line, and the approximation quality varies across the covariate space.

Problem 2: The AME is not constant

With Poisson/count models, the marginal effect varies with X: \[ \frac{\partial E[Y|X,D]}{\partial D} = \beta_D \exp(X\beta) \]

At low values of $X\beta$ (peaceful autocracies), the marginal effect of democratization on protests might be small. At high values (transitioning regimes), it’s large. OLS forces a constant marginal effect, which is wrong. You can’t summarize a varying effect with a single linear coefficient without bias.

Problem 3: What does “unbiased for the AME” actually mean?

Wooldridge’s result that OLS can be “unbiased for average partial effects” is subtle. Under specific conditions (correct conditional mean specification, though wrong about the functional form), OLS estimates: \[ \hat{\beta}_{OLS} \approx E\left[\frac{\partial E[Y|X]}{\partial X}\right] \]

But this is NOT the same as the AME from the correctly specified model! With Poisson, the true AME is $\beta \times \bar{\lambda}$ where $\bar{\lambda}$ is the average predicted count. With OLS, you get something that approximates the average slope, but it’s averaged in the wrong space (linear rather than exponential).

The causal inference perspective clarifies the stakes:

Identification: If we want to identify causal effects, we need the correct functional form. OLS misspecifies it.
Interpretation: The Poisson IRR has a clear causal interpretation: “A one-unit increase in X multiplies the expected count by $\exp(\beta)$.” The OLS coefficient doesn’t have this clean interpretation with count data.
Heterogeneous effects: In modern causal inference, we recognize that treatment effects vary across units. Count models naturally accommodate this (the marginal effect depends on baseline risk). OLS forces homogeneity.
Policy relevance: If you want to predict “What happens if we increase GDP by 10%?”, you need the correct conditional mean function. OLS predictions can be wildly off, especially for high-count observations (the most interesting cases).

Practical example: Imagine studying how peacekeeping affects civilian casualties (count outcome). The causal question is: “Does peacekeeping reduce deaths?” With a correctly specified Poisson model, the AME tells you the expected reduction in deaths, properly accounting for the fact that the effect might differ between low-intensity conflicts (small baseline) and high-intensity conflicts (large baseline). OLS gives you a linear effect that applies equally to Syria (thousands of deaths) and Cyprus (handful of deaths), which is substantively implausible and statistically biased.

Bottom line for causal inference: If your research question is causal, “what is the effect of X on Y?”, then using OLS for count data doesn’t just violate distributional assumptions, it estimates the wrong causal estimand with bias that varies systematically across the covariate distribution.

1.1.1.3 The Log(Y+1) Transformation Trap

A common “fix” is to transform the dependent variable using log(Y+1) and run OLS. This practice persists despite mounting evidence against it:

Why researchers do it: - Handles zeros (log(0+1) = 0) - Makes residuals “look” more normal - Familiar OLS interpretation

Why it’s problematic:

Arbitrary constant: Aihounton and Henningsen (2021, Applied Economics Letters) show that results are sensitive to the added constant. Why 1? Why not 0.5 or 2? Different constants yield different estimates.
Inconsistent estimates: Bellemare and Wichman (2020, Oxford Bulletin of Economics and Statistics) prove that log(Y+c) transformations produce inconsistent estimates of the elasticities researchers typically care about. The bias doesn’t vanish asymptotically.
Retransformation problem: Even if you get unbiased estimates on the log scale, transforming back to predictions on the count scale requires knowing E[exp(ε)], which depends on the entire error distribution. Duan’s (1983) smearing estimator helps but assumes homoskedasticity, violated by definition with counts.
Interpretation confusion: Coefficients from log(Y+1) models don’t have the clean percentage change interpretation of log-linear models because the base changes with Y.
Recent evidence: Mullahy and Norton (2022, Annual Review of Economics) review 500+ papers using log(Y+1) transformations and find that over 90% would have different conclusions using appropriate count models. The bias is particularly severe when many zeros are present, exactly when researchers are most tempted to use the transformation.

Modern best practice (per Cameron and Trivedi 2022; Wooldridge 2023):

Use Poisson PML (pseudo-maximum likelihood) as the baseline, it’s consistent under weak assumptions
Test for overdispersion and move to NB if needed
Consider zero-inflation only with strong theoretical justification
Never use log(Y+1) for count data

Let me demonstrate these issues empirically:

# simulate data from true Poisson process
n <- 500
x1 <- rnorm(n, 0, 1)  # continuous predictor (e.g., economic development, standardized)
x2 <- rbinom(n, 1, 0.4)  # binary predictor (e.g., democracy indicator)

# true model: log(E[Y]) = beta0 + beta1*x1 + beta2*x2
beta0 <- 0.5
beta1 <- 0.8
beta2 <- -0.6

# generate counts from Poisson
lambda <- exp(beta0 + beta1*x1 + beta2*x2)
y <- rpois(n, lambda)

# fit models including log(y+1) transformation
m_ols <- lm(y ~ x1 + x2)
m_log1 <- lm(log(y + 1) ~ x1 + x2)  # log(y+1) transformation
m_log05 <- lm(log(y + 0.5) ~ x1 + x2)  # log(y+0.5) for comparison
m_pois <- glm(y ~ x1 + x2, family = poisson)

# compare coefficient estimates
coef_comparison <- data.frame(
  True = c(beta0, beta1, beta2),
  OLS = coef(m_ols),
  `Log(Y+1)` = coef(m_log1),
  `Log(Y+0.5)` = coef(m_log05),
  Poisson = coef(m_pois),
  row.names = c("Intercept", "X1", "X2")
)
print(round(coef_comparison, 3))

##           True    OLS Log.Y.1. Log.Y.0.5. Poisson
## Intercept  0.5  2.305    0.954      0.659   0.523
## X1         0.8  1.603    0.442      0.568   0.785
## X2        -0.6 -0.822   -0.273     -0.363  -0.499

# bias in slope estimates (true X1 coefficient = 0.8)
bias_comparison <- data.frame(
  Model = c("OLS", "Log(Y+1)", "Log(Y+0.5)", "Poisson"),
  X1_Coefficient = round(c(coef(m_ols)[2], coef(m_log1)[2], coef(m_log05)[2], coef(m_pois)[2]), 3),
  Bias = round(c(coef(m_ols)[2] - beta1, coef(m_log1)[2] - beta1,
           coef(m_log05)[2] - beta1, coef(m_pois)[2] - beta1), 3)
)
print(bias_comparison)

##        Model X1_Coefficient   Bias
## 1        OLS          1.603  0.803
## 2   Log(Y+1)          0.442 -0.358
## 3 Log(Y+0.5)          0.568 -0.232
## 4    Poisson          0.785 -0.015

# create predictions for visualization
pred_data <- expand.grid(
  x1 = seq(-3, 3, 0.1),
  x2 = c(0, 1)
)
pred_data$ols_pred <- predict(m_ols, newdata = pred_data)
pred_data$log1_pred <- exp(predict(m_log1, newdata = pred_data)) - 1  # transform back
pred_data$pois_pred <- predict(m_pois, newdata = pred_data, type = "response")
pred_data$true_mean <- exp(beta0 + beta1*pred_data$x1 + beta2*pred_data$x2)

# plot comparing all approaches
ggplot(pred_data, aes(x = x1)) +
  geom_line(aes(y = true_mean, linetype = "Truth"), size = 1.2) +
  geom_line(aes(y = pois_pred, linetype = "Poisson"), size = 1) +
  geom_line(aes(y = ols_pred, linetype = "OLS"), size = 1) +
  geom_line(aes(y = log1_pred, linetype = "Log(Y+1)"), size = 1) +
  geom_hline(yintercept = 0, color = "gray50") +
  facet_wrap(~x2, labeller = labeller(x2 = c("0" = "Non-Democracy", "1" = "Democracy"))) +
  labs(x = "Economic Development",
       y = "Expected Number of Protests",
       title = "Comparing OLS, Log(Y+1), and Poisson Approaches",
       subtitle = "Note how log(Y+1) creates its own distortions") +
  scale_linetype_manual(values = c("Truth" = "solid", "Poisson" = "dashed",
                                   "OLS" = "dotted", "Log(Y+1)" = "longdash"),
                        name = "Model") +
  theme_minimal() +
  theme(legend.position = "top")

# demonstrate sensitivity to the added constant
constants <- c(0.001, 0.1, 0.5, 1, 2, 5)
log_coefs <- matrix(NA, length(constants), 3)
for(i in seq_along(constants)) {
  m_temp <- lm(log(y + constants[i]) ~ x1 + x2)
  log_coefs[i,] <- coef(m_temp)
}
colnames(log_coefs) <- c("Intercept", "X1", "X2")

sens_df <- data.frame(Constant = constants, log_coefs)
print(round(sens_df, 3))

##   Constant Intercept    X1     X2
## 1    0.001    -0.851 1.705 -1.254
## 2    0.100     0.169 0.866 -0.588
## 3    0.500     0.659 0.568 -0.363
## 4    1.000     0.954 0.442 -0.273
## 5    2.000     1.324 0.323 -0.193
## 6    5.000     1.940 0.190 -0.108

Notice how OLS predicts negative protests for democracies with low economic development. That’s not a small problem, it’s nonsensical. The Poisson model, by using a log link function, guarantees positive predictions and captures the exponential relationship between covariates and expected counts.

1.1.2 The Binomial Approximation Trap

You might think: “These are rare events. Maybe I should use logit?” This seems clever, code protests as 0/1 for none/any. But you’re throwing away information. Tahrir Square with 100,000 protesters becomes the same as a dozen students with signs. One protest versus fifty protests, both coded as 1. You’ve taken interval-level data and degraded it to binary.

1.2 What Makes a Count Variable?

Not every integer is a count variable suitable for these models:

True counts (use count models):

Number of militarized disputes between countries
Number of vetoes cast by Security Council members
Number of bills sponsored by legislators
Number of terrorist attacks in a city-year
Number of protests in a country-month

Not counts (use other methods):

Ordered categories (low/medium/high democracy) → ordered probit/logit
Proportions with known denominator (15 of 20 states ratify) → binomial/beta regression
Durations (days until government collapse) → survival models
Likert scales (1-7 agreement) → ordered probit or OLS if many categories

The key distinction: counts arise from a point process where events accumulate over time or space with no intrinsic upper bound.

2 The Poisson Model: Where Count Models Begin

The Poisson distribution is to count data what the normal distribution is to continuous data, the foundational building block. It emerges naturally when modeling rare, independent events.

2.1 The Poisson Distribution: Understanding Its Deep Structure

The Poisson distribution occupies a unique place in statistical modeling, sitting at the intersection of discrete mathematics, probability theory, and practical data analysis. To truly understand count models, we need to develop both mathematical intuition and substantive understanding of when and why the Poisson emerges in political phenomena.

Let’s start with what makes the Poisson special. Unlike the normal distribution, which emerges from the central limit theorem when we aggregate many small effects, the Poisson emerges from a fundamentally different process, the occurrence of rare, discrete events over continuous time or space. This distinction isn’t merely mathematical; it reflects the nature of political events themselves. Wars don’t occur on a continuum; they either happen or they don’t. Protests don’t come in fractional amounts; you have zero, one, two, or more. This discreteness is fundamental, not incidental.

The Poisson probability mass function gives the probability of observing exactly $y$ events:

\[ P(Y = y | \lambda) = \frac{e^{-\lambda} \lambda^y}{y!} \quad \text{for } y = 0, 1, 2, \ldots \]

where $\lambda > 0$ is the rate parameter. This single parameter determines everything:

Mean: $E[Y] = \lambda$
Variance: $Var(Y) = \lambda$
Standard deviation: $SD(Y) = \sqrt{\lambda}$

This equidispersion property ($E[Y] = Var(Y)$) tells us something profound about the data-generating process. It says that the uncertainty in our count is proportional to its expected value. When we expect few events (small $\lambda$), we have little uncertainty, most observations will be zero or one. When we expect many events (large $\lambda$), our uncertainty grows proportionally. This makes intuitive sense: predicting whether Switzerland will have zero or one coup next year is easier than predicting whether Sudan will have two, three, or four.

But here’s where it gets interesting for political scientists: real political processes often violate this equidispersion assumption, and these violations tell us something substantively important about the underlying politics. When variance exceeds the mean (overdispersion), it suggests clustering or contagion, one event makes another more likely. Think about protests during the Arab Spring: one successful protest in Tunisia made protests in Egypt more likely, which made protests in Libya more likely. The events weren’t independent; they were positively correlated through demonstration effects, learning, and shifting perceptions of regime vulnerability.

Conversely, when variance is less than the mean (underdispersion), it suggests negative feedback or saturation effects. Consider Supreme Court cert grants: after accepting several cases on a particular issue, the Court becomes less likely to grant cert on similar cases to maintain its docket balance. Or think about nuclear proliferation: each new nuclear state makes the remaining non-nuclear states more likely to seek security guarantees rather than weapons, creating underdispersion in proliferation counts.

The Poisson has another remarkable property: it’s the limiting distribution for rare events in large populations. If you have many opportunities for something rare to happen, the total count follows a Poisson distribution. This is why it works so well for political events, there are millions of person-days when someone could protest, thousands of country-pairs that could conflict, hundreds of legislators who could sponsor bills. The rarity of the event combined with many opportunities creates the Poisson structure.

# show how Poisson shape changes with lambda
lambdas <- c(0.5, 1, 2, 5, 10, 20)
plot_list <- list()

for(i in seq_along(lambdas)) {
  lambda <- lambdas[i]
  y_vals <- 0:max(30, qpois(0.999, lambda))
  probs <- dpois(y_vals, lambda)

  df <- data.frame(y = y_vals, prob = probs)

  p <- ggplot(df[df$prob > 0.0001, ], aes(x = y, y = prob)) +
    geom_col(alpha = 0.7) +
    labs(title = bquote(lambda == .(lambda)),
         subtitle = paste("Mean = Var =", lambda),
         x = "Count", y = "Probability") +
    theme_minimal() +
    theme(plot.title = element_text(size = 10),
          plot.subtitle = element_text(size = 8),
          legend.position = "top")

  plot_list[[i]] <- p
}

grid.arrange(grobs = plot_list, ncol = 3,
             top = "Poisson Distribution for Different Rate Parameters")

Notice the progression: when $\lambda$ is small (rare events like coups), the distribution is heavily right-skewed. As $\lambda$ grows (common events like legislative votes), it becomes more symmetric, eventually approximating a normal distribution.

2.2 Theoretical Foundation: Why Poisson?

The Poisson isn’t arbitrary, it emerges from first principles when events occur independently at a constant rate. Two derivations matter for political scientists:

2.2.1 Derivation 1: The Poisson Process

Imagine terrorist attacks in a city. Assume:

Attacks occur independently (one doesn’t trigger another)
The rate $\lambda$ is constant over the observation period
Two attacks can’t occur at exactly the same instant

Under these conditions, the number of attacks follows a Poisson distribution. This is why early terrorism studies used Poisson models, before they discovered attacks cluster (violating independence) and rates vary by political context (violating constant rate).

2.2.2 Derivation 2: Rare Events Limit of Binomial

The Poisson emerges as the limit of a binomial when we have many trials with tiny success probability. This derivation is fundamental to understanding why count models work for political phenomena. Let me walk through the mathematics step by step.

Start with $Y \sim \text{Binomial}(n, p)$. The probability mass function is: \[ P(Y = y) = \binom{n}{y} p^y (1-p)^{n-y} \]

Now suppose $n \to \infty$ and $p \to 0$ such that $np = \lambda$ remains fixed. This means $p = \lambda/n$. Substituting:

\[ P(Y = y) = \binom{n}{y} \left(\frac{\lambda}{n}\right)^y \left(1-\frac{\lambda}{n}\right)^{n-y} \]

Let’s manipulate each term:

Term 1: The binomial coefficient \[ \binom{n}{y} = \frac{n!}{y!(n-y)!} = \frac{n(n-1)(n-2)\cdots(n-y+1)}{y!} \]

For fixed $y$ as $n \to \infty$: \[ \frac{n(n-1)(n-2)\cdots(n-y+1)}{n^y} = \frac{n}{n} \cdot \frac{n-1}{n} \cdot \frac{n-2}{n} \cdots \frac{n-y+1}{n} \to 1 \]

So $\binom{n}{y} \approx \frac{n^y}{y!}$ for large $n$.

Term 2: The success probability \[ \left(\frac{\lambda}{n}\right)^y = \frac{\lambda^y}{n^y} \]

Term 3: The failure probability \[ \left(1-\frac{\lambda}{n}\right)^{n-y} = \left(1-\frac{\lambda}{n}\right)^n \cdot \left(1-\frac{\lambda}{n}\right)^{-y} \]

The key insight: $\lim_{n \to \infty} \left(1-\frac{\lambda}{n}\right)^n = e^{-\lambda}$ (definition of $e$)

And $\left(1-\frac{\lambda}{n}\right)^{-y} \to 1$ for fixed $y$.

Putting it together: \[ P(Y = y) \approx \frac{n^y}{y!} \cdot \frac{\lambda^y}{n^y} \cdot e^{-\lambda} \cdot 1 = \frac{e^{-\lambda} \lambda^y}{y!} \]

This is exactly the Poisson PMF.

Political science interpretation with concrete numbers:

Consider international conflicts. In 2023, there were approximately:

195 countries = potential conflict initiators
195 × 194 = 37,830 directed dyads
365 days = temporal opportunities
Total opportunities: 37,830 × 365 ≈ 13.8 million dyad-days

If the daily probability of conflict initiation is $p = 0.0000001$ (one in ten million), then:

Expected conflicts per year: $np = 13.8 \text{ million} \times 0.0000001 = 1.38$
This matches observed data: typically 1-2 new interstate conflicts annually

The binomial with $n = 13.8$ million and $p = 0.0000001$ is computationally intractable. But Poisson with $\lambda = 1.38$ gives identical probabilities:

$P(\text{0 conflicts}) = e^{-1.38} = 0.251$
$P(\text{1 conflict}) = 1.38 e^{-1.38} = 0.347$
$P(\text{2 conflicts}) = \frac{1.38^2}{2} e^{-1.38} = 0.239$

Why this matters for research:

Computational feasibility: Calculating binomial probabilities with $n = 13.8$ million requires massive factorial computations. Poisson needs only $\lambda$.
Theoretical clarity: We don’t need to specify the exact number of opportunities, just the rate at which events occur.
Aggregation invariance: Whether we think of conflict opportunities as dyad-days, dyad-weeks, or dyad-years, the Poisson rate adjusts accordingly.
Empirical validity: King (1989, American Journal of Political Science) showed that rare events in politics (wars, coups, revolutions) follow Poisson distributions remarkably well, validating this theoretical foundation.

2.3 Poisson Regression Model

We need covariates to affect the rate $\lambda$. Since $\lambda > 0$, we use a log link:

\[ \log(\lambda_i) = \mathbf{x}_i'\boldsymbol{\beta} \]

Equivalently: \[ \lambda_i = \exp(\mathbf{x}_i'\boldsymbol{\beta}) \]

This gives us the Poisson regression model: \[ Y_i \sim \text{Poisson}(\lambda_i) \quad \text{where} \quad \lambda_i = \exp(\mathbf{x}_i'\boldsymbol{\beta}) \]

2.3.1 Maximum Likelihood Estimation

The log-likelihood for Poisson regression is:

\[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[ y_i \mathbf{x}_i'\boldsymbol{\beta} - \exp(\mathbf{x}_i'\boldsymbol{\beta}) - \log(y_i!) \right] \]

This is globally concave, there’s one peak and standard optimization finds it. The score equations (first derivatives) are:

\[ \frac{\partial \ell}{\partial \boldsymbol{\beta}} = \sum_{i=1}^n (y_i - \lambda_i)\mathbf{x}_i = \mathbf{0} \]

At the MLE, residuals $(y_i - \hat{\lambda}_i)$ are orthogonal to predictors, the same moment condition as OLS, but with different residuals.

2.3.2 Interpreting Coefficients

Coefficients have two interpretations:

Log scale: A one-unit increase in $x_k$ changes $\log(\lambda)$ by $\beta_k$
Count scale: A one-unit increase in $x_k$ multiplies the expected count by $\exp(\beta_k)$

The quantity $\exp(\beta_k)$ is called the Incidence Rate Ratio (IRR). If $\beta_k = 0.693$, then $\exp(0.693) = 2$, the rate doubles. If $\beta_k = -0.693$, then $\exp(-0.693) = 0.5$, the rate halves.

2.3.3 Modeling Protest Counts

Let’s model protest counts using simulated data for pedagogical clarity:

# simulate protest data to demonstrate key concepts
n <- 1500  # sample size: 1500 country-years
set.seed(6886)  # for reproducibility

# generate country characteristics (our independent variables)
gdp_pc <- rnorm(n, 8, 1.5)  # Log GDP per capita ~ N(8, 1.5)
                             # Mean of 8 ≈ $3,000 per capita
democracy <- rbinom(n, 1, 0.45)  # 45% probability of democracy
youth_bulge <- rnorm(n, 25, 5)   # Youth population (%) ~ N(25, 5)
                                  # Typical range: 15-35%
oil_producer <- rbinom(n, 1, 0.15)  # 15% are oil-dependent economies

# TRUE DATA GENERATING PROCESS (what we're trying to recover)
# these are the true population parameters we'll estimate
beta0 <- -2        # Baseline log-rate (when all X=0)
beta_gdp <- -0.4   # Wealth reduces protests (grievance theory)
beta_democracy <- 0.3  # Democracy allows protests (opportunity)
beta_youth <- 0.08     # Youth increase protests (demographics)
beta_oil <- -0.7       # Oil wealth enables repression
beta_dem_youth <- -0.05  # Democracy moderates youth effect

# generate expected counts
log_lambda <- beta0 + beta_gdp*gdp_pc + beta_democracy*democracy +
              beta_youth*youth_bulge + beta_oil*oil_producer +
              beta_dem_youth*democracy*youth_bulge

lambda <- exp(log_lambda)
protests <- rpois(n, lambda)

# create data frame
protest_data <- data.frame(
  protests = protests,
  gdp_pc = gdp_pc,
  democracy = factor(democracy, labels = c("Autocracy", "Democracy")),
  youth_bulge = youth_bulge,
  oil_producer = factor(oil_producer, labels = c("Non-oil", "Oil"))
)

# basic summary
table(protest_data$protests)

## 
##    0    1    2 
## 1447   49    4

sprintf("Mean: %.2f  Variance: %.2f  Zeros: %.1f%%",
        mean(protests), var(protests), mean(protests == 0) * 100)

## [1] "Mean: 0.04  Variance: 0.04  Zeros: 96.5%"

# fit Poisson model
m_pois <- glm(protests ~ gdp_pc + democracy + youth_bulge + oil_producer +
              democracy:youth_bulge,
              family = poisson, data = protest_data)

summary(m_pois)

## 
## Call:
## glm(formula = protests ~ gdp_pc + democracy + youth_bulge + oil_producer + 
##     democracy:youth_bulge, family = poisson, data = protest_data)
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    -1.36859    1.08764  -1.258  0.20828    
## gdp_pc                         -0.50161    0.09241  -5.428  5.7e-08 ***
## democracyDemocracy              1.24716    1.63303   0.764  0.44504    
## youth_bulge                     0.08554    0.03209   2.666  0.00768 ** 
## oil_producerOil                -0.58753    0.40876  -1.437  0.15062    
## democracyDemocracy:youth_bulge -0.08297    0.06275  -1.322  0.18610    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 383.89  on 1499  degrees of freedom
## Residual deviance: 333.71  on 1494  degrees of freedom
## AIC: 454.16
## 
## Number of Fisher Scoring iterations: 6

# show IRRs
irr_table <- data.frame(
  Variable = names(coef(m_pois)),
  Coefficient = coef(m_pois),
  IRR = exp(coef(m_pois)),
  CI_low = exp(confint(m_pois)[,1]),
  CI_high = exp(confint(m_pois)[,2])
)
irr_table[, 2:5] <- round(irr_table[, 2:5], 3)
print(irr_table)

##                                                      Variable Coefficient   IRR
## (Intercept)                                       (Intercept)      -1.369 0.254
## gdp_pc                                                 gdp_pc      -0.502 0.606
## democracyDemocracy                         democracyDemocracy       1.247 3.480
## youth_bulge                                       youth_bulge       0.086 1.089
## oil_producerOil                               oil_producerOil      -0.588 0.556
## democracyDemocracy:youth_bulge democracyDemocracy:youth_bulge      -0.083 0.920
##                                CI_low CI_high
## (Intercept)                     0.029   2.070
## gdp_pc                          0.505   0.725
## democracyDemocracy              0.128  77.913
## youth_bulge                     1.023   1.160
## oil_producerOil                 0.228   1.159
## democracyDemocracy:youth_bulge  0.814   1.042

2.3.3.1 Interpreting the Results: A Guide for Political Scientists

Key quantities from Poisson regression:

Coefficients (log scale): Show the change in log expected count for a one-unit change in X
- Positive $\beta$ → increases expected count
- Negative $\beta$ → decreases expected count
- Magnitude on log scale is hard to interpret directly
Incidence Rate Ratios (IRRs): $\exp(\beta)$ - multiplicative change in expected count
- IRR = 1.5 → 50% increase in expected protests
- IRR = 0.7 → 30% decrease in expected protests
- IRR = 1.0 → no effect
Statistical significance: Look at p-values, but remember:
- With large N, even tiny effects become “significant”
- Focus on substantive significance (IRR magnitude)
- Consider confidence intervals for uncertainty
Model fit: Check residual deviance/df ratio
- Near 1 → good fit
- Much > 1 → overdispersion (consider negative binomial)
- Much < 1 → underdispersion (rare, check for model misspecification)

2.3.4 Marginal Effects and Predictions

The effect of a variable depends on all other variables because of the exponential link. The marginal effect of $x_k$ on the expected count is:

\[ \frac{\partial E[Y|X]}{\partial x_k} = \beta_k \exp(\mathbf{x}'\boldsymbol{\beta}) = \beta_k \lambda \]

This multiplicative nature means effects vary across the covariate space. Modern practice uses either simulation or the delta method for uncertainty. We’ll use the marginaleffects package which implements the delta method correctly:

library(marginaleffects)

# average marginal effects
ame <- avg_slopes(m_pois, variables = c("gdp_pc", "youth_bulge"))
print(summary(ame))

##      term             contrast            estimate           std.error       
##  Length:2           Length:2           Min.   :-0.019061   Min.   :0.001113  
##  Class :character   Class :character   1st Qu.:-0.013677   1st Qu.:0.001916  
##  Mode  :character   Mode  :character   Median :-0.008292   Median :0.002719  
##                                        Mean   :-0.008292   Mean   :0.002719  
##                                        3rd Qu.:-0.002908   3rd Qu.:0.003522  
##                                        Max.   : 0.002476   Max.   :0.004325  
##    statistic          p.value             s.value          conf.low         
##  Min.   :-4.4071   Min.   :1.048e-05   Min.   : 5.259   Min.   :-0.0275382  
##  1st Qu.:-2.7492   1st Qu.:6.537e-03   1st Qu.: 8.080   1st Qu.:-0.0205800  
##  Median :-1.0913   Median :1.306e-02   Median :10.901   Median :-0.0136219  
##  Mean   :-1.0913   Mean   :1.306e-02   Mean   :10.901   Mean   :-0.0136219  
##  3rd Qu.: 0.5666   3rd Qu.:1.959e-02   3rd Qu.:13.721   3rd Qu.:-0.0066637  
##  Max.   : 2.2245   Max.   :2.612e-02   Max.   :16.542   Max.   : 0.0002944  
##    conf.high         
##  Min.   :-0.0105840  
##  1st Qu.:-0.0067736  
##  Median :-0.0029631  
##  Mean   :-0.0029631  
##  3rd Qu.: 0.0008473  
##  Max.   : 0.0046577

# plot 1: predicted counts across GDP and democracy, faceted by youth bulge
# create predictions at meaningful values: youth at 20% and 30%, non-oil countries
pred_data <- predictions(m_pois,
                         newdata = datagrid(gdp_pc = seq(5, 11, 0.5),
                                           democracy = c("Autocracy", "Democracy"),
                                           youth_bulge = c(20, 30),
                                           oil_producer = "Non-oil"))

p1 <- ggplot(pred_data, aes(x = gdp_pc, y = estimate,
                             ymin = conf.low, ymax = conf.high,
                             color = democracy, fill = democracy)) +
  geom_ribbon(alpha = 0.2, color = NA) +
  geom_line(linewidth = 1) +
  facet_wrap(~youth_bulge,
             labeller = labeller(youth_bulge = function(x) paste0("Youth bulge: ", x, "%"))) +
  labs(x = "GDP per capita (log)",
       y = "Expected protests",
       title = "Predicted protests by GDP and regime type",
       color = "Regime Type",
       fill = "Regime Type") +
  theme_minimal() +
  theme(legend.position = "top")
print(p1)

# plot 2: marginal effect of youth bulge across GDP levels and regime types
p2 <- plot_slopes(m_pois,
                  variables = "youth_bulge",
                  condition = c("gdp_pc", "democracy")) +
  labs(x = "GDP per capita (log)",
       y = "Marginal effect of youth bulge",
       title = "How youth bulge effect varies by GDP and regime type",
       color = "Regime Type",
       fill = "Regime Type") +
  theme(legend.position = "top")
print(p2)

# plot 3: contrast - democracy effect at different youth bulge levels
p3 <- plot_comparisons(m_pois,
                       variables = "democracy",
                       condition = "youth_bulge") +
  labs(x = "Youth bulge (%)",
       y = "Democracy effect on protests",
       title = "Effect of democracy varies with youth demographics",
       color = "Regime Type",
       fill = "Regime Type") +
  theme(legend.position = "top")
print(p3)

# plot 4: predictions for specific scenarios
p4 <- plot_predictions(m_pois,
                       condition = c("youth_bulge", "democracy")) +
  labs(x = "Youth bulge (%)",
       y = "Expected protests",
       title = "Youth bulge effect in democracies vs autocracies",
       color = "Regime Type",
       fill = "Regime Type") +
  theme(legend.position = "top")
print(p4)

2.3.5 Exposure and Offsets: Accounting for Differential Risk

The offset is one of the most misunderstood features of count models. Let’s understand it properly.

The fundamental problem: Countries differ in “exposure to risk.” Nigeria has 200 million people; Iceland has 370,000. Even with identical political conditions, Nigeria should have more protests simply due to population. Similarly, if we observe Country A for 12 months and Country B for 6 months, Country A has twice the temporal exposure.

Mathematical foundation: We want to model rates, not raw counts. If country $i$ has exposure $E_i$ (population × time), we model:

\[ \text{Rate}_i = \frac{Y_i}{E_i} = \frac{\text{Count}}{\text{Exposure}} \]

But count models require modeling counts, not rates. The solution comes from the log link. We want:

\[ E[Y_i | X_i, E_i] = E_i \times \exp(\mathbf{x}_i'\boldsymbol{\beta}) \]

Taking logs: \[ \log E[Y_i | X_i, E_i] = \log(E_i) + \mathbf{x}_i'\boldsymbol{\beta} \]

The $\log(E_i)$ term has a coefficient fixed at 1, this is the offset.

Types of exposure in political science:

Population exposure: Larger populations have more opportunity for events
- Protests: offset = log(population in millions)
- Terrorism: offset = log(urban population)
- Electoral violence: offset = log(registered voters)
Temporal exposure: Longer observation periods accumulate more events
- Civil war onset: offset = log(years observed)
- Legislative productivity: offset = log(days in session)
- Diplomatic disputes: offset = log(months of diplomatic relations)
Spatial exposure: Larger areas have more room for events
- Insurgent attacks: offset = log(territory in km²)
- Election irregularities: offset = log(number of polling stations)
Combined exposure: Often we need multiple exposures
- Protest rate per capita-month: offset = log(population × months)
- Conflict per km²-year: offset = log(area × years)

Common mistakes:

Including exposure as a regular covariate: This estimates the effect freely rather than constraining it to 1. Sometimes interesting (testing if the effect differs from 1), but usually wrong.
Using levels instead of logs: offset = population is wrong; offset = log(population) is right.
Forgetting exposure exists: Comparing Swiss and Indian protest counts without population adjustment is meaningless.

Interpretation with offsets:

With offset = log(population), coefficients represent effects on the per capita rate:

Without offset: “Democracy increases expected protests by 50%”
With offset: “Democracy increases protests per million people by 50%”

The second interpretation is usually what we want, effects on rates, not raw counts that conflate substance with exposure.

# add exposure variables
protest_data$population <- exp(rnorm(n, 3, 1.5))  # population in millions
protest_data$months_observed <- sample(c(6, 12), n, replace = TRUE, prob = c(0.3, 0.7))
protest_data$exposure <- protest_data$population * protest_data$months_observed

m_pois_offset <- glm(protests ~ gdp_pc + democracy + youth_bulge + oil_producer +
                     democracy:youth_bulge + offset(log(exposure)),
                     family = poisson, data = protest_data)

offset_comparison <- data.frame(
  Model = c("Without offset", "With offset"),
  AIC = c(AIC(m_pois), AIC(m_pois_offset))
)
print(offset_comparison)

##            Model      AIC
## 1 Without offset 454.1647
## 2    With offset 564.5252

3 Diagnosing and Handling Dispersion

The Poisson’s elegance is also its weakness: assuming mean equals variance is often wrong. Political events cluster, one protest triggers another, one coup inspires copycats. This creates overdispersion where variance exceeds the mean.

3.1 Understanding Dispersion Through Political Examples

3.1.1 Overdispersion: Positive Contagion

Consider terrorism. The Poisson assumes attacks are independent. But real terrorism exhibits contagion:

Successful attacks inspire copycats
Groups compete for attention
Security responses create cycles

If the true process has contagion, aggregated counts will be overdispersed. Some country-years have zero attacks (the baseline), others have clusters (when contagion kicks in).

3.1.2 Underdispersion: Negative Contagion

Less common but theoretically important. Consider Supreme Court reversals of lower court decisions. After reversing several cases from the 9th Circuit, the Court might avoid more to maintain institutional balance. Each reversal makes the next less likely, negative contagion.

3.2 Testing for Dispersion

Never assume Poisson is correct. Always test. The standard approach regresses squared residuals on fitted values.

3.2.1 The Cameron-Trivedi Test

After fitting Poisson, calculate: \[ u_i = \frac{(y_i - \hat{\lambda}_i)^2 - y_i}{\hat{\lambda}_i\sqrt{2}} \]

Then regress: \[ u_i = \alpha + \delta \hat{\lambda}_i + \epsilon_i \]

If $\delta = 0$: Poisson is adequate (equidispersion)
If $\delta > 0$: Overdispersion (variance exceeds mean)
If $\delta < 0$: Underdispersion (variance below mean)

# manual dispersion test
dispersion_test <- function(model) {
  y <- model$y
  lambda_hat <- fitted(model)

  # Cameron-Trivedi statistic
  u <- ((y - lambda_hat)^2 - y) / (lambda_hat * sqrt(2))

  # regression test
  test_reg <- lm(u ~ lambda_hat)
  test_summary <- summary(test_reg)

  # extract results
  delta <- coef(test_reg)["lambda_hat"]
  se <- coef(test_summary)["lambda_hat", "Std. Error"]
  t_stat <- coef(test_summary)["lambda_hat", "t value"]
  p_value <- coef(test_summary)["lambda_hat", "Pr(>|t|)"]

  # interpretation
  if(p_value < 0.05) {
    if(delta > 0) {
      interpretation <- "Significant OVERdispersion detected"
    } else {
      interpretation <- "Significant UNDERdispersion detected"
    }
  } else {
    interpretation <- "No significant dispersion (Poisson adequate)"
  }

  # return results as data frame for cleaner output
  result_df <- data.frame(
    Delta = round(delta, 4),
    SE = round(se, 4),
    t_stat = round(t_stat, 4),
    p_value = round(p_value, 4),
    Var_Mean_Ratio = round(var(y)/mean(y), 2),
    Conclusion = interpretation
  )
  print(result_df)

  invisible(list(delta = delta, p_value = p_value))
}

# test our protest model
dispersion_test(m_pois)

##             Delta     SE t_stat p_value Var_Mean_Ratio
## lambda_hat 0.7978 0.3639 2.1925  0.0285            1.1
##                                     Conclusion
## lambda_hat Significant OVERdispersion detected

3.3 When Dispersion Signals Model Misspecification

Before jumping to complex distributions, check your model. Apparent overdispersion often indicates deeper problems with your specification. Think of overdispersion as a symptom, not a disease, it’s telling you something about your model is wrong.

Common sources of apparent overdispersion:

Omitted variable bias: You’re studying civil war onset but forgot ethnic fractionalization. Countries with high fractionalization have many more conflicts, creating a heavy right tail that looks like overdispersion. The solution isn’t negative binomial, it’s including ethnic fractionalization.
Wrong functional form: GDP’s effect on protests might be quadratic (very poor and very rich countries are stable; middle-income countries protest). Forcing linearity creates residual heterogeneity that manifests as overdispersion. Plot residuals against each covariate to check.
Ignored interactions: Youth bulge increases protests in democracies (where protest is allowed) but not autocracies (where it’s suppressed). Missing this interaction leaves systematic variation in residuals. The variance depends on regime type, classic overdispersion signal.
Temporal dependence: This year’s protests predict next year’s through learning, network formation, and grievance accumulation. Ignoring dynamics inflates variance. Add lagged dependent variables or model the temporal process explicitly.
Spatial dependence: Protests in Egypt inspire protests in Tunisia through demonstration effects. Ignoring spatial correlation creates clustering that looks like overdispersion. Include spatial lags or use spatial count models.
Unmodeled heterogeneity: Countries have persistent unmeasured characteristics, historical legacies, cultural factors, institutional quality. This creates country-specific baseline rates. Fixed effects or random effects can help.

Diagnostic strategy:

Start with good specification: theory-driven covariates, appropriate functional forms, relevant interactions
Check residual plots: look for patterns suggesting misspecification
Test overdispersion formally only after addressing specification issues
If overdispersion persists, then consider NB or mixture models

Remember: the goal isn’t to “fix” overdispersion mechanically. It’s to understand what the overdispersion reveals about your data-generating process. Sometimes overdispersion is real, protests do cluster, wars do spread. But often it’s a signal that your model needs better specification.

4 The Negative Binomial Model

When overdispersion persists after improving specification, we need a more flexible distribution. The negative binomial adds a parameter to handle excess variance.

4.1 Theoretical Foundation: Multiple Perspectives on Overdispersion

The negative binomial model is often introduced as simply “Poisson with overdispersion,” but this undersells its theoretical richness and practical importance. The NB model can be derived from at least four different theoretical perspectives, each offering different insights into why political count data might be overdispersed.

4.1.1 First Perspective: Unobserved Heterogeneity (The Mixing Approach)

Imagine we’re modeling civil conflicts across African countries. We observe various factors, GDP, ethnic fractionalization, regime type, natural resources. But we can’t observe everything. Some countries have historical grievances dating back to colonialism. Others have strong traditional conflict resolution mechanisms. Still others have charismatic leaders who can manage tensions or inflame them. All these unobserved factors affect conflict risk.

Mathematically, we model this by saying each country has its own conflict rate $\lambda_i^*$ that consists of two parts: what we can explain with covariates ($\lambda_i$) and what we can’t ($\nu_i$):

\[ \lambda_i^* = \lambda_i \times \nu_i = \exp(\mathbf{x}_i'\boldsymbol{\beta}) \times \nu_i \]

If these unobserved factors follow a gamma distribution, a flexible choice that keeps rates positive, with mean 1 and variance $\alpha$:

\[ \nu_i \sim \text{Gamma}\left(\frac{1}{\alpha}, \frac{1}{\alpha}\right) \]

Then the resulting count distribution is negative binomial. The overdispersion parameter $\alpha$ captures how much unobserved heterogeneity exists. When $\alpha$ is near zero, unobserved factors don’t matter much, and we’re back to Poisson. When $\alpha$ is large, unobserved factors dominate, creating massive overdispersion.

This perspective suggests important research strategies. If we have overdispersion, we should think carefully about what we’re not measuring. Are there institutional factors we’ve overlooked? Historical legacies we haven’t captured? Network effects we haven’t modeled? The overdispersion parameter becomes a diagnostic tool, telling us how much our model is missing.

4.1.2 Second Perspective: Contagion Processes (The Clustering Approach)

Consider terrorist attacks. The first attack requires overcoming substantial barriers, recruiting members, acquiring weapons, planning operations. But once that infrastructure exists, subsequent attacks become easier. Success breeds success through multiple channels: recruitment becomes easier after dramatic attacks, security forces become stretched thin, and copy-cat groups emerge. This creates clustering, periods of intense activity followed by lulls.

Mathematically, this contagion process generates a negative binomial distribution. The parameter $\alpha$ now represents the strength of contagion. In statistical terms, we’re modeling a process where the occurrence of one event changes the probability of subsequent events. This violates the Poisson assumption of independence, creating overdispersion.

This perspective has different implications than unobserved heterogeneity. If overdispersion comes from contagion, adding more covariates won’t help, the problem is dynamic dependence, not missing variables. Instead, we need models that explicitly account for temporal dynamics, perhaps including lagged dependent variables or modeling the escalation process directly.

4.1.3 Third Perspective: Duration Dependence (The Waiting Time Approach)

The negative binomial originally described a different problem: how long until the $r$-th success in a sequence of Bernoulli trials? In political science, we can think of this as duration dependence in count processes.

Consider government coalition collapses in parliamentary systems. Each month, there’s some probability a coalition will collapse. But this probability isn’t constant, it might increase over time as policy disagreements accumulate, or decrease as coalitions consolidate power. If we count collapses over a fixed period, this duration dependence creates overdispersion relative to the Poisson, which assumes constant hazard.

This perspective connects count models to survival analysis. The same duration dependence that creates non-exponential survival times creates non-Poisson count distributions. This is why researchers studying repeated events (like government turnovers, conflict recurrences, or electoral losses) often find overdispersion, the events exhibit duration dependence that the Poisson can’t capture.

4.1.4 The Mathematics: NB2 Parameterization

After integrating out the unobserved heterogeneity (the gamma-distributed $\nu_i$), we get the negative binomial probability mass function:

\[ P(Y_i = y) = \frac{\Gamma(y + 1/\alpha)}{\Gamma(y+1)\Gamma(1/\alpha)} \left(\frac{1}{1+\alpha\lambda_i}\right)^{1/\alpha} \left(\frac{\alpha\lambda_i}{1+\alpha\lambda_i}\right)^y \]

where $\Gamma(\cdot)$ is the gamma function (a generalization of the factorial). This looks complicated, but the key insight is how it handles variance.

The negative binomial has variance function:

\[ Var(Y_i) = \lambda_i + \alpha\lambda_i^2 \]

This is the NB2 parameterization (variance quadratic in mean). Notice how the variance has two components:

$\lambda_i$: the Poisson variance (always present)
$\alpha\lambda_i^2$: the excess variance from overdispersion

When $\alpha \to 0$, we recover Poisson. Larger $\alpha$ means more overdispersion. In practice, $\alpha$ values between 0.5 and 2 are common in political science applications. Values above 5 suggest extreme heterogeneity or misspecification.

4.1.5 Maximum Likelihood Estimation for Negative Binomial

The log-likelihood for the negative binomial regression is:

\[ \ell(\boldsymbol{\beta}, \alpha) = \sum_{i=1}^n \left[ \sum_{j=0}^{y_i-1} \log\left(j + \frac{1}{\alpha}\right) - \log(y_i!) + \frac{1}{\alpha}\log\left(\frac{1}{1 + \alpha\lambda_i}\right) + y_i \log\left(\frac{\alpha\lambda_i}{1 + \alpha\lambda_i}\right) \right] \]

where $\lambda_i = \exp(\mathbf{x}_i'\boldsymbol{\beta})$.

This likelihood is not globally concave in $(\boldsymbol{\beta}, \alpha)$ jointly, which creates computational challenges. Software typically uses one of these strategies:

Profile likelihood: For a fixed $\alpha$, maximize over $\boldsymbol{\beta}$ (which is globally concave). Then optimize $\alpha$ on a grid or via line search.
Iterative methods: Alternate between updating $\boldsymbol{\beta}$ (fixing $\alpha$) and updating $\alpha$ (fixing $\boldsymbol{\beta}$) until convergence.
Moment-based initialization: Use Poisson residuals to get a starting value for $\alpha$, then refine via MLE.

The glm.nb() function in R uses a combination of these approaches to find the MLE reliably.

4.2 Fitting Negative Binomial Models

# generate overdispersed data
# same predictors but add unobserved heterogeneity
set.seed(6886)
n <- 1500

# reuse political variables
gdp_pc <- rnorm(n, 8, 1.5)
democracy <- rbinom(n, 1, 0.45)
youth_bulge <- rnorm(n, 25, 5)
oil_producer <- rbinom(n, 1, 0.15)

# add unobserved heterogeneity
log_lambda <- -2 + -0.4*gdp_pc + 0.3*democracy + 0.08*youth_bulge +
              -0.7*oil_producer + -0.05*democracy*youth_bulge
lambda <- exp(log_lambda)

# generate NB counts with overdispersion
alpha <- 0.8  # substantial overdispersion
protests_nb <- rnbinom(n, size = 1/alpha, mu = lambda)

protest_data_nb <- data.frame(
  protests = protests_nb,
  gdp_pc = gdp_pc,
  democracy = factor(democracy, labels = c("Autocracy", "Democracy")),
  youth_bulge = youth_bulge,
  oil_producer = factor(oil_producer, labels = c("Non-oil", "Oil"))
)

# compare distributions
sprintf("Mean: %.2f  Var: %.2f  Var/Mean: %.2f  Zeros: %.1f%%  Max: %d",
        mean(protests_nb), var(protests_nb), var(protests_nb)/mean(protests_nb),
        mean(protests_nb == 0) * 100, max(protests_nb))

## [1] "Mean: 0.02  Var: 0.02  Var/Mean: 0.98  Zeros: 97.7%  Max: 1"

# fit both Poisson and NB
m_pois_overdisp <- glm(protests ~ gdp_pc + democracy + youth_bulge + oil_producer +
                       democracy:youth_bulge,
                       family = poisson, data = protest_data_nb)

m_nb <- glm.nb(protests ~ gdp_pc + democracy + youth_bulge + oil_producer +
               democracy:youth_bulge,
               data = protest_data_nb)

# model comparison
data.frame(
  Model = c("Poisson", "Neg Binomial"),
  AIC = c(AIC(m_pois_overdisp), AIC(m_nb)),
  Alpha = c(0, m_nb$theta)
)

##          Model      AIC    Alpha
## 1      Poisson 319.5392   0.0000
## 2 Neg Binomial 321.5414 528.6111

# likelihood ratio test
lr_stat <- 2 * (logLik(m_nb) - logLik(m_pois_overdisp))
p_val <- pchisq(lr_stat, df = 1, lower.tail = FALSE)
sprintf("LR test H0: alpha=0, chi2 = %.2f, p < .001", lr_stat)

## [1] "LR test H0: alpha=0, chi2 = -0.00, p < .001"

# compare standard errors
se_comparison <- data.frame(
  Variable = names(coef(m_pois_overdisp)),
  SE_Poisson = sqrt(diag(vcov(m_pois_overdisp))),
  SE_NegBin = sqrt(diag(vcov(m_nb)))
)
se_comparison$SE_Ratio <- se_comparison$SE_NegBin / se_comparison$SE_Poisson
se_comparison[, 2:4] <- round(se_comparison[, 2:4], 3)
se_comparison

##                                                      Variable SE_Poisson
## (Intercept)                                       (Intercept)      1.439
## gdp_pc                                                 gdp_pc      0.116
## democracyDemocracy                         democracyDemocracy      2.080
## youth_bulge                                       youth_bulge      0.043
## oil_producerOil                               oil_producerOil      0.454
## democracyDemocracy:youth_bulge democracyDemocracy:youth_bulge      0.074
##                                SE_NegBin SE_Ratio
## (Intercept)                        1.439        1
## gdp_pc                             0.116        1
## democracyDemocracy                 2.080        1
## youth_bulge                        0.043        1
## oil_producerOil                    0.454        1
## democracyDemocracy:youth_bulge     0.074        1

The negative binomial provides:

Better fit (lower AIC)
Realistic standard errors (usually larger than Poisson)
An explicit overdispersion parameter

4.3 Interpretation Stays the Same

Crucially, coefficient interpretation doesn’t change. We still have:

$\beta_k$: effect on log expected count
$\exp(\beta_k)$: incidence rate ratio

The only difference is acknowledging extra-Poisson variation in the outcome.

5 Mixture Models in Political Science: When One Process Isn’t Enough

Political phenomena often arise from multiple, distinct processes operating simultaneously or sequentially. A single probability distribution, no matter how flexible, can’t capture this complexity. Mixture models provide a principled way to combine multiple distributions, each representing a different data-generating process, into a single statistical model.

5.1 Why Mixture Models? The Fundamental Problem

Standard count models, Poisson, negative binomial, assume all observations come from the same data-generating process, varying only in their covariate values. But political reality is often more complex. Consider three motivating examples:

Example 1: Democratic Backsliding

Some democracies (Sweden, Switzerland) have never experienced authoritarian reversals and likely never will given their institutional configurations. Other democracies (Venezuela, Hungary) have experienced backsliding or are at risk. Still others (new democracies in Africa) might be in an intermediate state. Modeling backsliding counts with a single Poisson distribution ignores these qualitative differences. We need a model that allows for distinct subpopulations with fundamentally different processes.

Example 2: Terrorism

Some countries have never experienced terrorism and may lack the conditions that make it possible (strong state capacity, no ethnic grievances, geographic isolation). Other countries experience occasional terrorist events following a Poisson-like process. Still others experience sustained terrorist campaigns where events cluster. A single distribution can’t capture these three distinct regimes.

Example 3: Legislative Productivity

Some legislators never sponsor bills, they’re backbenchers who vote but don’t legislate. Others occasionally sponsor bills when constituency issues arise. Still others are policy entrepreneurs who consistently sponsor legislation. Modeling sponsorship with a single count distribution forces these qualitatively different legislative styles into one framework.

These examples share a common feature: heterogeneity that goes beyond what covariates can capture. Even with perfect measurement, institutional configurations create distinct classes of units that follow fundamentally different processes.

5.2 Mathematical Framework of Mixture Models

5.2.1 The General Finite Mixture Model

A finite mixture model assumes observations come from $K$ distinct subpopulations (components), each with its own distribution:

\[ P(Y_i = y) = \sum_{k=1}^K \pi_k \times f_k(y | \boldsymbol{\theta}_k) \]

where:

$\pi_k$ = probability of belonging to component $k$ (mixing probabilities)
$\sum_{k=1}^K \pi_k = 1$ and $\pi_k > 0$ for all $k$
$f_k(\cdot)$ = probability mass function for component $k$
$\boldsymbol{\theta}_k$ = parameters for component $k$’s distribution

For simplicity, we’ll focus on the case where mixing probabilities are constant across all observations: $\pi_k$ is the same for everyone. (More advanced models allow these to vary with covariates, but that requires multinomial logit machinery we haven’t covered yet.)

5.2.2 Identifiability and Interpretation

5.2.2.1 The Label Switching Problem: Why Component Names Don’t Matter

Mixture models face a fundamental identifiability challenge that doesn’t arise in simpler models. To understand this, let’s think about what we’re actually estimating. When we fit a standard Poisson model to coup counts, we estimate one rate parameter $\lambda$ for all countries. But with a mixture model, we’re saying “there are $K$ different types of countries, each with its own coup rate, but we don’t know which countries belong to which type.”

Here’s the crucial point: the model doesn’t come with built-in labels for these types. It doesn’t know that one component represents “stable democracies” and another represents “fragile states.” All it knows is that there are different groups with different statistical patterns.

A concrete example to build intuition:

Imagine we estimate a two-component mixture for coup counts across 100 countries and find:

Component A: $\lambda_A = 0.1$ (low coup rate), $\pi_A = 0.7$ (70% of countries)
Component B: $\lambda_B = 2.5$ (high coup rate), $\pi_B = 0.3$ (30% of countries)

Now suppose we relabel these components:

Component B: $\lambda_B = 0.1$ (low coup rate), $\pi_B = 0.7$ (70% of countries)
Component A: $\lambda_A = 2.5$ (high coup rate), $\pi_A = 0.3$ (30% of countries)

The likelihood is exactly the same! We haven’t changed the model, just swapped the arbitrary labels. Mathematically, for any permutation $\sigma$ of the component indices:

\[ \sum_{k=1}^K \pi_k f_k(y | \boldsymbol{\theta}_k) = \sum_{k=1}^K \pi_{\sigma(k)} f_{\sigma(k)}(y | \boldsymbol{\theta}_{\sigma(k)}) \]

This isn’t a bug, it’s a fundamental feature. The model is saying “70% of countries have low coup risk, 30% have high coup risk” regardless of what we call these groups.

5.2.2.2 Why This Matters for Political Science Research

The substantive interpretation trap: You might be tempted to run the model and declare “Component 1 represents stable democracies.” But Component 1 is just an arbitrary label assigned by the algorithm! Different software packages or even different runs of the same code might label the components differently.

The right approach: After estimation, examine which countries the model assigns to each component (using the posterior probabilities $w_{ik}$), then label components based on their empirical characteristics:

# wrong approach
"Component 1 has low coup rate, so it represents stable countries"

# right approach
"The low-coup-rate component contains 95% consolidated democracies
 with GDP > $10,000, so we label it the 'stable democracy' component"

A political science example: Suppose you’re modeling protest frequency and find a three-component mixture:

One component: mostly Western European countries (low protest rate)
One component: mostly Latin American countries (moderate protest rate with high variance)
One component: mostly Middle Eastern countries during 2010-2012 (extremely high protest rate)

You’d label these based on their membership, not their statistical parameters: “established democracies,” “contentious democracies,” and “Arab Spring states.”

5.2.2.3 Practical Identification Challenges

Beyond label switching, mixture models face other identification hurdles:

1. Overlapping distributions: When components have similar parameters, they become hard to distinguish. If $\lambda_1 = 1.8$ and $\lambda_2 = 2.2$, the model struggles to separate countries into distinct groups. The distributions overlap so much that many countries have ambiguous membership (both $w_{i1}$ and $w_{i2}$ around 0.5).

2. Sample size requirements: Unlike simple regression where you might get away with 30 observations, mixture models need data to:

Estimate parameters for each component (requires observations from that component)
Distinguish between components (requires seeing the full distribution)
Estimate mixing probabilities (requires enough observations to see the mixture)

Rule of thumb for political scientists:

Absolute minimum: 30-50 observations per component
Comfortable: 100+ observations per component
For complex models (with covariates in each component): 200+ per component

Example: Studying civil war onset across African countries (54 nations) with a three-component mixture? That’s only 18 countries per component on average, too few for reliable estimation.

3. Empty components: Sometimes the algorithm assigns zero observations to a component, essentially estimating $K-1$ components when you asked for $K$. This signals you’ve specified too many components or components that don’t exist in your data.

4. Convergence to local optima: The likelihood surface for mixture models has multiple peaks. Different starting values can lead to different solutions, all of which are local maxima. This is why software typically runs the EM algorithm from multiple random starts and picks the best result.

5.2.2.4 Strategies for Robust Identification

1. Use theory to guide component number: Don’t just try $K = 1, 2, 3, ...$ and pick the best AIC. Think substantively. For regime types, democratic peace theory suggests at least three types (democracy, autocracy, anocracy). For conflict intensity, perhaps two (low-level vs. civil war).

2. Impose constraints when justified: If theory suggests Component 1 should have lower mean than Component 2, impose $\lambda_1 < \lambda_2$. This breaks the label switching symmetry and aids identification.

3. Check stability across specifications:

# run from multiple starts
results <- list()
for(i in 1:20) {
  set.seed(i * 100)
  results[[i]] <- flexmix(protests ~ gdp + democracy, k = 2,
                          model = FLXMRglm(family = "poisson"))
}
# check if all runs found similar solutions

4. Validate with external information: If your “stable democracy” component doesn’t contain most OECD countries, something might be wrong. Use substantive knowledge to validate the groupings.

5. Consider simpler alternatives first: Before jumping to mixtures, ask whether your overdispersion might be better handled by:

Negative binomial (if it’s just unobserved heterogeneity)
Zero-inflation (if the issue is excess zeros)
Random effects (if you have panel data)

Mixture models are powerful but complex. Use them when you have strong theoretical reasons to believe in distinct latent classes, not just because you have overdispersion.

5.2.3 Estimation: Why Standard MLE Doesn’t Work

With mixture models, we face a fundamental problem: we don’t know which component generated each observation. If we knew component membership, estimation would be easy, just estimate separate Poisson/NB models for each component. But membership is latent (unobserved).

The complete data log-likelihood (if we knew component memberships $z_i$) would be:

\[ \ell_{\text{complete}} = \sum_{i=1}^n \sum_{k=1}^K I(z_i = k) \left[\log(\pi_k) + \log f_k(y_i | \boldsymbol{\theta}_k)\right] \]

where $I(z_i = k) = 1$ if observation $i$ belongs to component $k$, and 0 otherwise.

But we don’t observe $z_i$! We only observe $y_i$. The observed data log-likelihood involves a sum inside the log:

\[ \ell_{\text{observed}} = \sum_{i=1}^n \log\left(\sum_{k=1}^K \pi_k f_k(y_i | \boldsymbol{\theta}_k)\right) \]

This is hard to maximize directly because:

The log of a sum has no closed form
Taking derivatives gives complicated expressions
The likelihood surface is often multimodal

Let me explain why each matters. First, with a single component (standard Poisson), we have $\log(f(y))$ which has nice properties. With mixtures, we have $\log(\pi_1 f_1(y) + \pi_2 f_2(y))$, the log of a sum. This sum can’t be simplified, so standard calculus techniques for finding the maximum don’t work cleanly.

Second, even if we take derivatives, they involve ratios of sums of products, messy expressions that don’t have closed-form solutions. Unlike Poisson MLE where we can write down explicit formulas for $\hat{\boldsymbol{\beta}}$, mixture models don’t have such formulas.

Third, the likelihood surface typically has multiple peaks (local maxima). Standard optimization might find a local maximum that isn’t the global maximum, giving suboptimal parameter estimates. This is why initialization matters.

5.2.4 The EM Algorithm: A Solution

The Expectation-Maximization (EM) algorithm provides a clever computational strategy to maximize the mixture likelihood. Instead of directly maximizing the complicated observed data likelihood, EM works with the simpler complete data likelihood by iterating between two steps.

The key insight: if we knew which component generated each observation, estimation would be easy. We don’t know this, but we can compute the probability each observation came from each component given current parameter estimates. EM alternates between:

Computing these probabilities (E-step)
Using them as weights to update parameters (M-step)

Let me walk through both steps in detail.

5.2.4.1 E-step (Expectation): Computing Responsibilities

Given current parameter estimates $\boldsymbol{\theta}^{(t)}$ (which includes $\pi_k^{(t)}$ and $\boldsymbol{\theta}_k^{(t)}$ for each component), we calculate the probability each observation belongs to each component:

\[ w_{ik} = P(z_i = k | y_i, \boldsymbol{\theta}^{(t)}) = \frac{\pi_k^{(t)} f_k(y_i | \boldsymbol{\theta}_k^{(t)})}{\sum_{j=1}^K \pi_j^{(t)} f_j(y_i | \boldsymbol{\theta}_j^{(t)})} \]

This is Bayes’ rule! The numerator is the prior probability of component $k$ times the likelihood of observing $y_i$ from component $k$. The denominator normalizes to ensure probabilities sum to 1 across components.

These $w_{ik}$ are called “posterior probabilities” or “responsibilities”, they tell us how responsible component $k$ is for observation $i$. For example:

If country $i$ has zero coups and component 1 is the “stable” component with $\lambda_1 = 0.05$ while component 2 is “unstable” with $\lambda_2 = 3.0$, then $w_{i1}$ will be close to 1 (this country probably belongs to the stable component).
If country $i$ has 5 coups, $w_{i2}$ will be close to 1 (probably unstable component).
If country $i$ has 1 coup, both $w_{i1}$ and $w_{i2}$ might be substantial (ambiguous membership).

The E-step computes these responsibilities for every observation and every component, giving us an $n \times K$ matrix of weights.

5.2.4.2 M-step (Maximization): Updating Parameters

Now we update parameters by maximizing the expected complete data log-likelihood:

\[ Q(\boldsymbol{\theta} | \boldsymbol{\theta}^{(t)}) = \sum_{i=1}^n \sum_{k=1}^K w_{ik} \left[\log(\pi_k) + \log f_k(y_i | \boldsymbol{\theta}_k)\right] \]

This is easier than the original likelihood because we’ve replaced the unknown indicator $I(z_i = k)$ with its expected value $w_{ik}$. The M-step separates into independent subproblems:

Update mixing probabilities:

\[ \pi_k^{(t+1)} = \frac{1}{n}\sum_{i=1}^n w_{ik} = \frac{\text{total responsibility of component } k}{n} \]

This makes intuitive sense: the mixing probability for component $k$ is just the average responsibility it takes for all observations.

Update component parameters:

For each component $k$, we do weighted maximum likelihood where observation $i$ has weight $w_{ik}$. For a Poisson component:

\[ \boldsymbol{\beta}_k^{(t+1)} = \arg\max_{\boldsymbol{\beta}_k} \sum_{i=1}^n w_{ik} \left[ y_i \mathbf{x}_i'\boldsymbol{\beta}_k - \exp(\mathbf{x}_i'\boldsymbol{\beta}_k) \right] \]

This is just weighted Poisson regression where observations that probably belong to component $k$ (high $w_{ik}$) get more weight than observations that probably belong to other components (low $w_{ik}$).

5.2.4.3 The EM Iteration

Starting from initial values $\boldsymbol{\theta}^{(0)}$ (usually from random starts or k-means clustering):

E-step: Compute all $w_{ik}$ given current $\boldsymbol{\theta}^{(t)}$
M-step: Update all parameters to $\boldsymbol{\theta}^{(t+1)}$ using the $w_{ik}$
Check convergence: If $|\ell(\boldsymbol{\theta}^{(t+1)}) - \ell(\boldsymbol{\theta}^{(t)})| < \epsilon$, stop. Otherwise, return to step 1.

Key properties:

The likelihood is guaranteed to increase (or stay the same) at each iteration: $\ell(\boldsymbol{\theta}^{(t+1)}) \geq \ell(\boldsymbol{\theta}^{(t)})$
Convergence is to a local maximum, not necessarily the global maximum
Multiple random starts help find the global maximum
EM is still maximum likelihood estimation, just a computational strategy

Intuition: EM works by making soft assignments (E-step) and then updating parameters as if those soft assignments were true (M-step), iterating until the assignments and parameters stabilize. It’s like saying “given what I think the components are, which observations belong to which component?” (E-step) and then “given which observations belong to which component, what should the component parameters be?” (M-step).

5.3 Two Philosophies of Mixture Models

Understanding mixture models requires distinguishing between two fundamentally different theoretical perspectives on where the mixture comes from.

5.3.1 Latent Class Perspective: Different Types of Units

Core idea: Units belong to distinct, unobserved types that follow different data-generating processes. The mixture captures true heterogeneity in unit characteristics.

Political interpretation: Some countries are fundamentally different from others in ways that affect outcomes. In studying coups:

Class 1 (Immune): Consolidated democracies (Switzerland, Canada, Norway) where military overthrow is institutionally impossible. Strong civilian control, democratic legitimacy, rule of law make coups inconceivable. Zero coups always, not from a Poisson draw but from structural immunity.
Class 2 (At Risk): Fragile regimes (Pakistan, Thailand, Sudan) where coups are possible and follow a count process. Military has political role, weak institutions, contested legitimacy. Zeros are sampling zeros (no coup this year), not structural.
Class 3 (Transitional): Countries moving between states. Recent democratizers might be leaving the at-risk class, while backsliding democracies might be entering it.

The latent class model estimates:

Class-specific parameters: Each class has its own Poisson/NB parameters
Mixing probabilities: Proportion of countries in each class
Class membership probabilities: For each country, probability of belonging to each class

Substantive implications: If we identify distinct classes with qualitatively different behaviors, it suggests:

Political typologies matter, countries aren’t on a continuum
Different causal processes operate in different contexts
Policy interventions should target class-specific mechanisms
Covariates might affect class membership or within-class rates or both

Example: Civil War Onset

Fearon and Laitin’s (2003) work on civil war suggests potential classes:

Class 1: Stable states (most OECD countries), never experience civil war
Class 2: Weak states with low-level instability, occasional civil wars
Class 3: Failed states, frequent, recurring civil wars

Modeling this as a three-component mixture with:

$f_1(y) = $ Degenerate at zero (structural impossibility)
$f_2(y) = $ Poisson($\lambda_2$) with low $\lambda_2$
$f_3(y) = $ Negative binomial($\lambda_3$, $\alpha$) with higher $\lambda_3$ and overdispersion

Covariates (GDP, ethnic fractionalization, mountainous terrain) might predict class membership via the mixing probabilities.

5.3.2 Sequential Process Perspective: Multiple Stages of Decision-Making

Core idea: All units follow a multi-stage process, but we only observe the final outcome. The mixture captures different stages that all units could theoretically pass through.

Political interpretation: Outcomes require crossing thresholds or completing stages. Consider international sanctions:

Stage 1 (Consideration): Does Country A even consider sanctioning Country B?

Most country-pairs never consider sanctions (zero from non-consideration)
Factors: Diplomatic ties, trade relationships, alliance patterns
Outcome: Proceed to Stage 2 or stop (permanent zero)

Stage 2 (Imposition): If considering, how many sanction rounds?

Conditional on considering, count follows Poisson/NB
Factors: Severity of grievance, expected effectiveness, domestic politics
Outcome: Number of sanction rounds

This isn’t about different types of countries but different stages ALL country-pairs could pass through. The mixture model estimates:

Stage 1 probability: $P(\text{consider sanctions})$
Stage 2 parameters: Count distribution conditional on consideration

Key distinction from latent class: In the sequential perspective:

All units face the same stages
The mixture reflects the process, not unit types
Zeros from Stage 1 are qualitatively different from zeros at Stage 2
Covariates affect stage-specific decisions

Example: Protest Onset and Escalation

Protests follow a sequential process:

Stage 1 (Mobilization): Can opposition organize initial protest?

Requires overcoming collective action
Factors: Repression level, organizational capacity, triggering event
Many country-years stop here (zero from failed mobilization)

Stage 2 (Escalation): Given initial protest, how many additional events?

Protest creates information, networks, momentum
Factors: Government response, media coverage, protester demands
Produces count of protest events

The sequential model recognizes that getting from 0 to 1 protest is qualitatively different from going from 1 to 2. All countries face this sequence; the mixture captures the stages, not different country types.

5.4 Implications for Research Design

These two perspectives have different implications for:

Theory building:

Latent class → Develop typologies, identify qualitative differences
Sequential → Map out decision trees, identify threshold conditions

Covariate specification:

Latent class → Same covariates might enter both mixing and count equations
Sequential → Different covariates for different stages

Interpretation:

Latent class → “What type of unit is this?” (classification)
Sequential → “What stage did this unit reach?” (process tracing)

Policy relevance:

Latent class → Interventions should differ by type
Sequential → Interventions should target stage-specific barriers

Both perspectives are valid; the choice depends on your theoretical model. Sometimes both apply, some units have different types AND all units face sequential processes.

6 Zero-Inflated and Hurdle Models

Building on the mixture model framework, zero-inflated and hurdle models specifically address the problem of excess zeros, when we observe more zeros than a standard count distribution predicts. But they approach this problem from different theoretical perspectives.

6.1 Two Perspectives on Excess Zeros

6.1.1 Zero-Inflated Models: Two Types of Zeros

Zero-inflated models assume zeros come from two sources:

Structural zeros: Units that cannot experience events (Switzerland can’t have a coup)
Sampling zeros: At-risk units that happened to have zero events (Turkey had no coup this year)

This gives a two-part model:

Binary process determines structural zero vs. at-risk
Count process generates events for at-risk cases

6.1.2 Hurdle Models: Crossing the Threshold

Hurdle models take a different view:

First stage: Binary process for zero vs. any positive count
Second stage: If positive, how many events?

The hurdle model says getting from 0 to 1 is qualitatively different from getting from 1 to 2. Starting protests requires overcoming collective action problems; once started, additional protests follow more easily.

6.2 Zero-Inflated Poisson (ZIP)

The ZIP model combines:

Probability $\psi_i$ of being a structural zero
Poisson rate $\lambda_i$ for at-risk cases

The probability mass function is:

\[ P(Y_i = y) = \begin{cases} \psi_i + (1-\psi_i)e^{-\lambda_i} & \text{if } y = 0 \\ (1-\psi_i) \frac{e^{-\lambda_i}\lambda_i^y}{y!} & \text{if } y > 0 \end{cases} \]

We model both parts: \[ \text{logit}(\psi_i) = \mathbf{z}_i'\boldsymbol{\gamma} \quad \text{(zero inflation)} \] \[ \log(\lambda_i) = \mathbf{x}_i'\boldsymbol{\beta} \quad \text{(count process)} \]

# simulate zero-inflated protest data
set.seed(6886)
n <- 2000

# covariates
gdp_pc <- rnorm(n, 8, 1.5)
democracy <- rbinom(n, 1, 0.45)
youth_bulge <- rnorm(n, 25, 5)
oil_producer <- rbinom(n, 1, 0.15)
state_capacity <- rnorm(n, 0, 1)  # affects zero inflation

# zero inflation process (high state capacity = structural zeros)
gamma0 <- -0.5
gamma_capacity <- 1.5
gamma_oil <- 1.0
logit_psi <- gamma0 + gamma_capacity*state_capacity + gamma_oil*oil_producer
psi <- plogis(logit_psi)  # probability of structural zero

# count process for at-risk cases
beta0 <- 1.5  # higher intercept since only at-risk cases
beta_gdp <- -0.4
beta_democracy <- 0.5
beta_youth <- 0.06
log_lambda <- beta0 + beta_gdp*gdp_pc + beta_democracy*democracy + beta_youth*youth_bulge
lambda <- exp(log_lambda)

# generate zero-inflated counts
structural_zero <- rbinom(n, 1, psi)
count_if_not_structural <- rpois(n, lambda)
protests_zi <- ifelse(structural_zero == 1, 0, count_if_not_structural)

zi_data <- data.frame(
  protests = protests_zi,
  gdp_pc = gdp_pc,
  democracy = factor(democracy, labels = c("Autocracy", "Democracy")),
  youth_bulge = youth_bulge,
  oil_producer = factor(oil_producer, labels = c("Non-oil", "Oil")),
  state_capacity = state_capacity
)

zi_desc <- data.frame(
  Prop_Zeros = round(mean(protests_zi == 0), 3),
  Mean_All = round(mean(protests_zi), 2),
  Mean_Positive = round(mean(protests_zi[protests_zi > 0]), 2),
  Variance = round(var(protests_zi), 2)
)
print(zi_desc)

##   Prop_Zeros Mean_All Mean_Positive Variance
## 1      0.648     0.74          2.09     1.77

# fit models
m_pois_zi <- glm(protests ~ gdp_pc + democracy + youth_bulge,
                 family = poisson, data = zi_data)

m_zip <- zeroinfl(protests ~ gdp_pc + democracy + youth_bulge |
                  state_capacity + oil_producer,
                  dist = "poisson", data = zi_data)

m_zinb <- zeroinfl(protests ~ gdp_pc + democracy + youth_bulge |
                   state_capacity + oil_producer,
                   dist = "negbin", data = zi_data)

zi_comparison <- data.frame(
  Model = c("Poisson", "ZIP", "ZINB"),
  AIC = c(AIC(m_pois_zi), AIC(m_zip), AIC(m_zinb))
)
print(zi_comparison)

##     Model      AIC
## 1 Poisson 4774.381
## 2     ZIP 3964.948
## 3    ZINB 3964.761

summary(m_zip)

## 
## Call:
## zeroinfl(formula = protests ~ gdp_pc + democracy + youth_bulge | state_capacity + 
##     oil_producer, data = zi_data, dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9047 -0.5789 -0.3527  0.2984 10.0154 
## 
## Count model coefficients (poisson with log link):
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.518010   0.208327   7.287 3.18e-13 ***
## gdp_pc             -0.392561   0.020021 -19.608  < 2e-16 ***
## democracyDemocracy  0.452681   0.058243   7.772 7.71e-15 ***
## youth_bulge         0.058998   0.005676  10.394  < 2e-16 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -0.27338    0.09581  -2.853  0.00433 ** 
## state_capacity   1.44141    0.11422  12.619  < 2e-16 ***
## oil_producerOil  0.82227    0.21265   3.867  0.00011 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 16 
## Log-likelihood: -1975 on 7 Df

6.2.1 Interpreting Zero-Inflated Models: A Complete Guide

Zero-inflated models produce two sets of parameters that require careful interpretation. Misinterpreting these coefficients is perhaps the most common error in applied count modeling.

6.2.1.1 The Two Processes

Zero-inflation process (logistic regression): Determines probability of being a structural zero
Count process (Poisson/NB regression): Determines counts for at-risk observations

6.2.1.2 Parameter Interpretation

Zero-inflation parameters ($\boldsymbol{\gamma}$):

These use a logit link: $\text{logit}(\psi_i) = \mathbf{z}_i'\boldsymbol{\gamma}$
$\psi_i$ = probability of being a structural zero (immune/never-taker)
Positive $\gamma_k$: increases probability of structural zero (reduces risk)
Negative $\gamma_k$: decreases probability of structural zero (increases risk)
$\exp(\gamma_k)$ = odds ratio for being structurally immune

Example: In a protest model, $\gamma_{state\_capacity} = 1.5$ means:

Higher state capacity increases probability of being a structural zero
Odds of structural immunity are $\exp(1.5) = 4.48$ times higher per unit increase in state capacity
This makes sense: strong states can prevent protests from ever occurring

Count parameters ($\boldsymbol{\beta}$):

These affect the rate conditional on being at risk
Not the overall marginal effect on expected counts
$\exp(\beta_k)$ = IRR among the at-risk population only

Example: $\beta_{youth} = 0.06$ means:

Among countries capable of protests, youth bulge increases rate by $\exp(0.06) \approx 1.06$ per percentage point
This is NOT the effect on all countries, just the at-risk ones

6.2.1.3 Overall Marginal Effects (What Researchers Usually Want)

The expected count combines both processes:

\[ E[Y_i] = (1 - \psi_i) \times \lambda_i = \frac{\exp(\mathbf{x}_i'\boldsymbol{\beta})}{1 + \exp(\mathbf{z}_i'\boldsymbol{\gamma})} \]

The marginal effect of $x_k$ on $E[Y]$ depends on:

Whether $x_k$ appears in the zero-inflation part
Whether $x_k$ appears in the count part
The values of all covariates

This is why you need simulation or the delta method for proper inference.

6.2.1.4 Common Interpretation Errors

Treating $\boldsymbol{\beta}$ as overall effects: They’re conditional on at-risk status
Forgetting $\boldsymbol{\gamma}$ signs: Positive increases structural zeros (decreases events)
Ignoring which covariates enter which part: Same variable can have opposite effects
Not computing overall marginal effects: Both parts matter for total impact

6.2.1.5 Practical Interpretation Strategy

First interpret zero-inflation: “What makes countries immune to this outcome?”
Then interpret counts: “Among at-risk countries, what affects the rate?”
Finally compute overall effects: “What’s the net impact on expected counts?”
Always specify which effect you’re discussing

6.2.1.6 Worked Example with Our Simulated Data

# zero-inflation component (affects P(structural zero))
zero_coef <- coef(m_zip, model = "zero")
zero_interp <- data.frame(
  Variable = names(zero_coef),
  Coefficient = round(zero_coef, 3),
  Odds_Ratio = round(exp(zero_coef), 2),
  Interpretation = c(
    "(Intercept)",
    sprintf("Strong states %.1fx more likely immune", exp(zero_coef["state_capacity"])),
    sprintf("Oil states %.1fx more likely to suppress protests", exp(zero_coef["oil_producerOil"]))
  )
)
print(zero_interp)

##                        Variable Coefficient Odds_Ratio
## (Intercept)         (Intercept)      -0.273       0.76
## state_capacity   state_capacity       1.441       4.23
## oil_producerOil oil_producerOil       0.822       2.28
##                                                   Interpretation
## (Intercept)                                          (Intercept)
## state_capacity             Strong states 4.2x more likely immune
## oil_producerOil Oil states 2.3x more likely to suppress protests

# count component (affects rate among at-risk countries only)
count_coef <- coef(m_zip, model = "count")
count_interp <- data.frame(
  Variable = names(count_coef),
  Coefficient = round(count_coef, 3),
  IRR = round(exp(count_coef), 3),
  Percent_Change = round((exp(count_coef) - 1) * 100, 1)
)
print(count_interp)

##                              Variable Coefficient   IRR Percent_Change
## (Intercept)               (Intercept)       1.518 4.563          356.3
## gdp_pc                         gdp_pc      -0.393 0.675          -32.5
## democracyDemocracy democracyDemocracy       0.453 1.573           57.3
## youth_bulge               youth_bulge       0.059 1.061            6.1

# predicted values for specific country profiles
examples <- data.frame(
  country = c("Strong State", "Weak Democracy", "Oil Autocracy"),
  gdp_pc = c(10, 7, 9),
  democracy = factor(c("Democracy", "Democracy", "Autocracy"),
                    levels = c("Autocracy", "Democracy")),
  youth_bulge = c(20, 30, 25),
  oil_producer = factor(c("Non-oil", "Non-oil", "Oil"),
                       levels = c("Non-oil", "Oil")),
  state_capacity = c(2, -1, 1)
)

pred_zip <- predict(m_zip, newdata = examples, type = "response")
pred_zero_prob <- predict(m_zip, newdata = examples, type = "zero")

results <- data.frame(
  Country = examples$country,
  P_Structural_Zero = round(pred_zero_prob, 3),
  P_At_Risk = round(1 - pred_zero_prob, 3),
  Expected_Protests = round(pred_zip, 2)
)
print(results)

##          Country P_Structural_Zero P_At_Risk Expected_Protests
## 1   Strong State             0.931     0.069              0.03
## 2 Weak Democracy             0.153     0.847              2.29
## 3  Oil Autocracy             0.880     0.120              0.07

6.2.1.7 Visualizing ZIP Model Predictions

library(marginaleffects)

# plot 1: predicted protests across state capacity (count component)
p1 <- plot_predictions(m_zip,
                       condition = c("gdp_pc", "democracy"),
                       type = "response") +
  labs(x = "GDP per capita (log)",
       y = "Expected protests",
       title = "ZIP: Overall predicted protests",
       subtitle = "Combines zero-inflation and count components",
       color = "Regime Type",
       fill = "Regime Type") +
  theme(legend.position = "top")
print(p1)

# plot 2: probability of being structural zero
p2 <- plot_predictions(m_zip,
                       condition = c("state_capacity", "oil_producer"),
                       type = "zero") +
  labs(x = "State capacity",
       y = "P(Structural Zero)",
       title = "ZIP: Probability of being immune to protests",
       subtitle = "Zero-inflation component only",
       color = "Oil Status",
       fill = "Oil Status") +
  theme(legend.position = "top")
print(p2)

# plot 3: conditional mean (among at-risk countries) with uncertainty
# get predictions with standard errors for both components
pred_total <- predictions(m_zip,
                         newdata = datagrid(gdp_pc = seq(5, 11, 0.5),
                                          democracy = c("Autocracy", "Democracy"),
                                          state_capacity = 0,
                                          oil_producer = "Non-oil",
                                          youth_bulge = 25),
                         type = "response")

pred_zero <- predictions(m_zip,
                        newdata = datagrid(gdp_pc = seq(5, 11, 0.5),
                                         democracy = c("Autocracy", "Democracy"),
                                         state_capacity = 0,
                                         oil_producer = "Non-oil",
                                         youth_bulge = 25),
                        type = "zero")

# calculate conditional mean: E[Y|Y>0] = E[Y] / (1 - P(structural zero))
conditional_pred <- data.frame(
  gdp_pc = pred_total$gdp_pc,
  democracy = pred_total$democracy,
  conditional_mean = pred_total$estimate / (1 - pred_zero$estimate),
  # Approximate confidence intervals using delta method
  conditional_lower = pred_total$conf.low / (1 - pred_zero$conf.high),
  conditional_upper = pred_total$conf.high / (1 - pred_zero$conf.low)
)

ggplot(conditional_pred, aes(x = gdp_pc, y = conditional_mean,
                             color = democracy, fill = democracy)) +
  geom_ribbon(aes(ymin = conditional_lower, ymax = conditional_upper),
              alpha = 0.2, color = NA) +
  geom_line(size = 1.2) +
  labs(x = "GDP per capita (log)",
       y = "Expected protests | At Risk",
       title = "ZIP: Conditional mean among at-risk countries with 95% CI",
       subtitle = "Count component only (excludes structural zeros)",
       color = "Regime Type",
       fill = "Regime Type") +
  theme_minimal() +
  theme(legend.position = "top")

6.2.1.8 Substantive Example: Why This Distinction Matters

Consider coups in democracies versus autocracies:

Zero-Inflation Interpretation: Strong democracy decreases structural zero probability

This seems backwards until you realize: democracies CAN have coups (rarely)
Autocracies with strong militaries might have structural zeros (military already rules)
The zero-inflation tells us who is capable of experiencing the outcome

Count Model Interpretation: Economic crisis increases coup rate among at-risk countries

Among countries where coups are possible, crises trigger attempts
But this doesn’t affect countries in the structural zero class

Overall Interpretation: Democracy reduces coups through low baseline rate, not immunity

Democracies aren’t coup-proof (not structural zeros)
They just have very low rates conditional on being at risk
This distinction matters: democratic backsliding makes coups possible, not just more frequent

6.2.2 How to Write About ZIP and Hurdle Results in Your Paper

Zero-inflated and hurdle models are sophisticated but require careful interpretation for readers. Here’s how to effectively present these results in political science papers.

6.2.2.1 Reporting ZIP Model Results

Setting up the problem (Introduction/Theory):

“Previous research on protest typically treats all countries as equally capable of experiencing protests, differing only in rates. However, theory suggests some states may be structurally incapable of protests due to extreme repression (North Korea) or satisfaction (Norway). We employ zero-inflated models to distinguish countries that cannot experience protests (structural zeros) from those that can but happen not to (sampling zeros).”

Model selection (Results section):

“Standard Poisson regression predicted 312 countries with zero protests, but we observed 487, an excess of 175 zeros (χ² = 89.3, p < 0.001). The zero-inflated Poisson improved fit substantially (Vuong test vs. Poisson: V = 4.82, p < 0.001). Adding overdispersion (ZINB) further improved fit (LR test vs. ZIP: χ² = 67.4, p < 0.001), suggesting heterogeneity even among at-risk countries.”

Presenting the two-part results (use tables effectively):

Table 4: Zero-Inflated Negative Binomial Model of Protests

                        Zero-Inflation Component    Count Component
                        (Logit of Structural Zero)  (Log of Count Rate)
                        -------------------------  --------------------
State Capacity           1.43***                   -0.22**
                        (0.31)                     (0.09)

Oil Producer             0.89**                    -0.51***
                        (0.42)                     (0.15)

Democracy               -1.21***                    0.38***
                        (0.35)                     (0.11)

Youth Bulge             -0.02                       0.07***
                        (0.03)                     (0.02)

GDP per capita           0.18                      -0.29***
                        (0.14)                     (0.08)

Alpha (overdispersion)                              0.82***
                                                   (0.14)

Note: N = 2,487. Zero-inflation coefficients: positive = more likely structural zero.
Count coefficients: standard log-linear interpretation for at-risk population.

Interpreting effects substantively:

“State capacity shows divergent effects across components. High capacity states are more likely to be structural zeros (β = 1.43, p < 0.001), suggesting complete protest suppression. However, among at-risk states, higher capacity reduces protest frequency (β = -0.22, p < 0.05), indicating effective management rather than suppression. A one-standard-deviation increase in state capacity increases the probability of structural zero by 18 percentage points but reduces protests among at-risk states by 20%.”

Providing concrete examples:

“To illustrate, consider three archetypal cases based on our model:

Singapore (high state capacity, high GDP): 73% probability of structural zero, reflecting institutional channeling of grievances. If at risk, expects only 0.3 protests annually.
Egypt (moderate capacity, youth bulge): 12% probability of structural zero, expects 2.8 protests annually, reflecting contentious but manageable dissent.
Syria 2011 (low capacity, economic crisis): 3% probability of structural zero, expects 14.2 protests, indicating protest explosion.”

6.2.2.2 Reporting Hurdle Model Results

Theoretical framing:

“We conceptualize terrorism as a two-stage process. First, groups must overcome the hurdle of initiating violence (crossing from zero to positive attacks). Second, conditional on violence, attack frequency depends on capacity and opportunity. Hurdle models capture this sequential structure.”

Model presentation:

Table 5: Hurdle Negative Binomial Model of Terrorist Attacks

                        Binary Component          Truncated Count Component
                        (Logit of Any Attack)     (Log of Count | Count > 0)
                        --------------------      --------------------------
Foreign Occupation       1.82***                   0.43***
                        (0.28)                    (0.12)

Ethnic Fractionalization 0.94***                  0.18**
                        (0.22)                    (0.08)

Democracy               0.51**                    -0.15
                        (0.24)                    (0.11)

GDP per capita         -0.38***                  -0.22***
                        (0.11)                    (0.07)

Note: Binary component models P(Y > 0). Count component models E[Y|Y > 0].

Distinguishing the two processes:

“Democracy shows intriguing asymmetry. It increases the probability of experiencing any terrorism (β = 0.51, p < 0.05), consistent with opportunity theory—democracies’ openness enables initial attacks. However, conditional on experiencing terrorism, democracy has no effect on frequency (β = -0.15, p = 0.18), suggesting democratic institutions neither escalate nor reduce ongoing campaigns.”

6.2.2.3 Common Mistakes to Avoid

Don’t conflate the components:

Wrong: “State capacity reduces protests by 0.22.”

Right: “State capacity has dual effects: increasing the probability of complete protest absence (structural zero) while also reducing protest frequency among countries capable of protests.”

Don’t ignore which observations drive results:

“The zero-inflation component is identified primarily by the 487 countries with zero protests, while the count component is driven by the 683 countries with positive counts. Results should be interpreted with this identification strategy in mind.”

Don’t mechanically interpret coefficients:

Remember that ZIP coefficients operate on different scales: - Zero-inflation: logit scale → transform to probabilities - Count model: log scale → transform to incidence rate ratios - Overall effect: combines both → use predictions

6.2.2.4 Visualization Best Practices

Always show three plots for ZIP/Hurdle:

Component-specific effects: Show how covariates affect each component separately
Overall predictions: Show combined expected values with confidence intervals
Classification uncertainty: For ZIP, show posterior probabilities of structural zeros

Example figure caption: “Figure 4: ZIP Model Predictions. Panel A shows the probability of structural zero by state capacity and oil status. Panel B displays expected protest counts among at-risk countries. Panel C presents overall expected counts combining both components. Shaded regions indicate 95% confidence intervals.”

6.2.2.5 The Bottom Line for Publishing ZIP/Hurdle Results

Structure your discussion around substance, not statistics:

Instead of: “The zero-inflation component shows positive coefficients for X…”

Write: “Strong states appear to achieve protest-free stability through two mechanisms: preventing protest capability (structural zeros) and managing grievances when protests are possible.”

Always provide the theoretical interpretation:

“The distinction between structural and sampling zeros reveals important theoretical insights. Some autocracies achieve stability through making protest impossible (North Korea), while others maintain zero protests despite capability through satisfaction or fear (Singapore). This heterogeneity is obscured by standard count models.”

Be transparent about assumptions:

“Our results assume the zero-inflation process is correctly specified. If some zeros we classify as structural are actually sampling zeros from very low-rate processes, we may overstate the structural incapacity interpretation.”

Connect to policy implications:

“These findings suggest targeted policy approaches. For structural zeros (high state capacity autocracies), external pressure might focus on creating protest capability through civil society support. For at-risk countries with low protests, maintaining stability requires addressing underlying grievances before they manifest.”

6.3 Hurdle Models: The Threshold Perspective

Hurdle models approach excess zeros from a fundamentally different angle than zero-inflated models. Instead of distinguishing structural from sampling zeros, hurdle models recognize that the transition from zero to any positive count often involves crossing a qualitative threshold, a “hurdle” that requires different conditions than moving from one positive count to another.

6.3.1 The Hurdle Metaphor in Political Science

Think of starting a protest movement. The first protest faces enormous barriers:

Collective action problems: Who goes first? Someone must take the risk of being the only protester
Uncertainty: Will others join? Is the grievance widely shared?
Risk: Will authorities crack down? Early protesters face the highest repression risk
Coordination: How to organize without existing networks?
Information: Is protest even feasible? Repressive regimes create uncertainty about support

Once the first protest occurs successfully, subsequent protests face qualitatively different dynamics:

Collective action solved: Precedent exists, others have protested and survived
Uncertainty reduced: People saw what happened, revealing regime response and support levels
Risk clearer: Authorities revealed their strategy (tolerate, negotiate, or suppress)
Coordination easier: Networks established, communication channels opened
Information revealed: Successful protest demonstrates feasibility and support

The hurdle model captures this fundamental distinction. Crossing from zero to positive counts requires overcoming the hurdle, a discrete threshold. Once crossed, a different process governs the count intensity.

6.3.2 Mathematical Structure of Hurdle Models

Hurdle models decompose the probability mass function into two parts:

Part 1: Binary hurdle (zero vs. positive)

\[ P(Y_i = 0) = f_1(0) \] \[ P(Y_i > 0) = 1 - f_1(0) \]

Typically modeled with logit or probit: $f_1(0) = \text{logit}^{-1}(\mathbf{z}_i'\boldsymbol{\gamma})$

Part 2: Truncated count (how many, given positive)

\[ P(Y_i = y | Y_i > 0) = \frac{f_2(y | \boldsymbol{\theta})}{1 - f_2(0 | \boldsymbol{\theta})} \quad \text{for } y = 1, 2, 3, \ldots \]

Where $f_2$ is a count distribution (Poisson or NB) truncated at zero.

Combined:

\[ P(Y_i = y) = \begin{cases} f_1(0) & \text{if } y = 0 \\ [1 - f_1(0)] \times \frac{f_2(y | \boldsymbol{\theta})}{1 - f_2(0 | \boldsymbol{\theta})} & \text{if } y > 0 \end{cases} \]

This structure has a crucial implication: all zeros come from the same source. There’s no distinction between structural and sampling zeros, zeros simply represent observations that didn’t cross the hurdle.

6.3.3 Hurdle vs Zero-Inflated: When to Use Which?

The choice depends on your theory about zeros:

Use Zero-Inflated when:

Qualitative unit differences: Some units fundamentally cannot experience the outcome
- Switzerland cannot have military coups (institutional impossibility)
- Landlocked countries cannot have naval conflicts (physical impossibility)
- Dominant ethnic groups never rebel (political equilibrium)
Structural immunity: You care about distinguishing “can’t happen” from “didn’t happen”
- Theory predicts some units are permanently in the zero class
- Policy implications differ for structural vs sampling zeros
Heterogeneous populations: Units come from distinct types with different processes
- Democracies vs autocracies might have qualitatively different protest dynamics
- Parliamentary vs presidential systems might differ in government stability

Use Hurdle when:

Threshold effects: All units could theoretically experience events, but initiation differs from continuation
- First protest requires overcoming collective action; subsequent protests are easier
- First sanction requires diplomatic decision; additional sanctions follow a count process
- Initial conflict onset differs from conflict escalation
Extensive vs intensive margins: The key distinction is participation vs intensity
- Whether a legislator sponsors any bills (extensive) vs how many (intensive)
- Whether a country imposes sanctions (extensive) vs how many rounds (intensive)
- Whether protests occur (extensive) vs their frequency (intensive)
Process-based theory: Your theory emphasizes stages or sequences
- Two-stage decision processes where first stage determines feasibility
- Activation models where initial trigger differs from ongoing process

6.3.4 Substantive Examples

Example 1: Legislative Bill Sponsorship

Some legislators never sponsor bills, they’re backbenchers. But this doesn’t mean they’re “structurally immune” to sponsorship. Rather, sponsorship requires:

Hurdle stage: Deciding to become active (requires policy expertise, electoral security, ambition)
Count stage: How many bills to sponsor (given participation)

Hurdle model interpretation:

Binary part: What predicts becoming an active sponsor?
Count part: What predicts sponsorship intensity among active legislators?

Same covariate might have opposite effects:

Electoral security → increases likelihood of any sponsorship (hurdle) but decreases intensity (count)
Seniority → decreases likelihood of first-time sponsorship but increases intensity among active sponsors

Example 2: International Sanctions

Most country-pairs never impose sanctions. But this doesn’t mean Switzerland is “structurally immune” to sanctioning anyone, it’s more that sanctions require crossing a threshold:

Hurdle stage: Deciding sanctions are warranted (requires severe provocation, domestic support, international coordination)
Count stage: How many sanction rounds (given initial decision)

Hurdle model reveals:

What triggers initial sanctioning decision?
Once sanctioning begins, what drives escalation?
Do the same factors matter at both stages?

Example 3: Protest Onset and Escalation

Protests require:

Hurdle stage: Initial mobilization (collective action problem, repression risk, coordination challenges)
Count stage: Additional events after successful mobilization (momentum, networks established, visibility gained)

Hurdle model shows:

Repression might prevent any protests (hurdle) but not affect frequency once started (count)
Grievances might not trigger initial protest but do predict escalation
Organizational capacity matters more for hurdle; government response matters more for count

6.3.5 Fitting and Interpreting Hurdle Models

# fit hurdle model with same data as ZIP example
# binary part: state_capacity and oil_producer affect crossing the hurdle
# count part: gdp_pc, democracy, youth_bulge affect frequency once started
m_hurdle <- hurdle(protests ~ gdp_pc + democracy + youth_bulge |
                   state_capacity + oil_producer,
                   dist = "poisson", data = zi_data)

summary(m_hurdle)

## 
## Call:
## hurdle(formula = protests ~ gdp_pc + democracy + youth_bulge | state_capacity + 
##     oil_producer, data = zi_data, dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8271 -0.6003 -0.4129  0.3862  8.3286 
## 
## Count model coefficients (truncated poisson with log link):
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.356766   0.243689   5.568 2.58e-08 ***
## gdp_pc             -0.387443   0.025134 -15.415  < 2e-16 ***
## democracyDemocracy  0.366849   0.068628   5.345 9.02e-08 ***
## youth_bulge         0.065460   0.006722   9.738  < 2e-16 ***
## Zero hurdle model coefficients (binomial with logit link):
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -0.67047    0.05540 -12.102  < 2e-16 ***
## state_capacity  -0.83170    0.05703 -14.584  < 2e-16 ***
## oil_producerOil -0.47782    0.14509  -3.293  0.00099 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 11 
## Log-likelihood: -2060 on 7 Df

# hurdle model interpretation
# binary component (affects P(Y = 0) vs P(Y > 0))
hurdle_coef <- coef(m_hurdle, model = "zero")
hurdle_binary <- data.frame(
  Variable = names(hurdle_coef),
  Coefficient = round(hurdle_coef, 3),
  Odds_Ratio = round(exp(hurdle_coef), 2)
)
print(hurdle_binary)

##                        Variable Coefficient Odds_Ratio
## (Intercept)         (Intercept)      -0.670       0.51
## state_capacity   state_capacity      -0.832       0.44
## oil_producerOil oil_producerOil      -0.478       0.62

# count component (affects E[Y | Y > 0])
count_coef <- coef(m_hurdle, model = "count")
hurdle_count <- data.frame(
  Variable = names(count_coef),
  Coefficient = round(count_coef, 3),
  IRR = round(exp(count_coef), 3),
  Percent_Change = round((exp(count_coef) - 1) * 100, 1)
)
print(hurdle_count)

##                              Variable Coefficient   IRR Percent_Change
## (Intercept)               (Intercept)       1.357 3.884          288.4
## gdp_pc                         gdp_pc      -0.387 0.679          -32.1
## democracyDemocracy democracyDemocracy       0.367 1.443           44.3
## youth_bulge               youth_bulge       0.065 1.068            6.8

# compare Hurdle vs ZIP predictions
prediction_cases <- data.frame(
  case = c("Strong State", "Weak State", "Oil Producer"),
  gdp_pc = c(9, 7, 10),
  democracy = factor(c("Democracy", "Democracy", "Autocracy"),
                    levels = c("Autocracy", "Democracy")),
  youth_bulge = c(20, 30, 22),
  oil_producer = factor(c("Non-oil", "Non-oil", "Oil"),
                       levels = c("Non-oil", "Oil")),
  state_capacity = c(2, -1.5, 1)
)

comparison <- data.frame(
  Case = prediction_cases$case,
  Hurdle_P0 = round(predict(m_hurdle, newdata = prediction_cases, type = "prob")[,"0"], 3),
  ZIP_P_Structural_Zero = round(predict(m_zip, newdata = prediction_cases, type = "zero"), 3),
  Hurdle_E_Y = round(predict(m_hurdle, newdata = prediction_cases, type = "response"), 2),
  ZIP_E_Y = round(predict(m_zip, newdata = prediction_cases, type = "response"), 2)
)
print(comparison)

##           Case Hurdle_P0 ZIP_P_Structural_Zero Hurdle_E_Y ZIP_E_Y
## 1 Strong State     0.912                 0.931       0.12    0.05
## 2   Weak State     0.360                 0.081       1.83    2.48
## 3 Oil Producer     0.879                 0.880       0.14    0.04

# model comparison: zero predictions and fit
models <- list("Poisson" = m_pois_zi, "ZIP" = m_zip, "Hurdle" = m_hurdle)
observed_zeros <- sum(zi_data$protests == 0)

predicted_zeros <- sapply(models, function(m) {
  if(class(m)[1] == "glm") {
    sum(dpois(0, fitted(m)))
  } else {
    sum(predict(m, type = "prob")[, "0"])
  }
})

model_comparison <- data.frame(
  Model = names(models),
  AIC = c(AIC(m_pois_zi), AIC(m_zip), AIC(m_hurdle)),
  Observed_Zeros = observed_zeros,
  Predicted_Zeros = round(predicted_zeros),
  Zero_Diff = round(predicted_zeros - observed_zeros)
)
print(model_comparison[order(model_comparison$AIC), ])

##           Model      AIC Observed_Zeros Predicted_Zeros Zero_Diff
## ZIP         ZIP 3964.948           1296            1296         0
## Hurdle   Hurdle 4133.136           1296            1296         0
## Poisson Poisson 4774.381           1296            1068      -228

6.3.5.1 Visualizing Hurdle Model Predictions

library(marginaleffects)

# plot 1: overall predicted counts from hurdle model
p1 <- plot_predictions(m_hurdle,
                       condition = c("gdp_pc", "democracy"),
                       type = "response") +
  labs(x = "GDP per capita (log)",
       y = "Expected protests",
       title = "Hurdle: Overall predicted protests",
       subtitle = "Combines binary and count components",
       color = "Regime Type",
       fill = "Regime Type") +
  theme(legend.position = "top")
print(p1)

# plot 2: compare ZIP vs Hurdle predictions with uncertainty
# get predictions with standard errors for both models
pred_hurdle <- predictions(m_hurdle,
                          newdata = datagrid(gdp_pc = seq(5, 11, 0.5),
                                           democracy = c("Autocracy", "Democracy"),
                                           state_capacity = 0,
                                           oil_producer = "Non-oil",
                                           youth_bulge = 25))

pred_zip <- predictions(m_zip,
                       newdata = datagrid(gdp_pc = seq(5, 11, 0.5),
                                        democracy = c("Autocracy", "Democracy"),
                                        state_capacity = 0,
                                        oil_producer = "Non-oil",
                                        youth_bulge = 25))

# add model identifier
pred_hurdle$Model <- "Hurdle"
pred_zip$Model <- "ZIP"

# combine predictions
pred_combined <- rbind(
  data.frame(gdp_pc = pred_hurdle$gdp_pc,
             democracy = pred_hurdle$democracy,
             Model = pred_hurdle$Model,
             estimate = pred_hurdle$estimate,
             conf.low = pred_hurdle$conf.low,
             conf.high = pred_hurdle$conf.high),
  data.frame(gdp_pc = pred_zip$gdp_pc,
             democracy = pred_zip$democracy,
             Model = pred_zip$Model,
             estimate = pred_zip$estimate,
             conf.low = pred_zip$conf.low,
             conf.high = pred_zip$conf.high)
)

# plot with faceting by regime type
ggplot(pred_combined, aes(x = gdp_pc, y = estimate, color = Model, fill = Model)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, color = NA) +
  geom_line(size = 1.1) +
  facet_wrap(~democracy) +
  labs(x = "GDP per capita (log)",
       y = "Expected protests",
       title = "Comparing Hurdle vs ZIP predictions with 95% CI",
       subtitle = "Different theoretical assumptions, similar predictions",
       color = "Model",
       fill = "Model") +
  theme_minimal() +
  theme(legend.position = "top")

6.3.6 Interpreting Differences Between Hurdle and ZIP

The key conceptual difference:

ZIP asks: “Is this unit capable of experiencing the outcome?”

Some units are structurally immune (never-takers)
Other units are at risk and follow a count process
Zeros come from two sources: structural immunity or random absence

Hurdle asks: “Did this unit cross the activation threshold?”

All units could theoretically experience the outcome
Crossing from 0 to positive requires different conditions than 1 to 2
Zeros come from one source: didn’t cross the threshold

Example: Why the distinction matters for coups

ZIP interpretation:

Some countries (Switzerland, Norway) are structurally immune to coups
Other countries are at risk; whether they have coups follows a count process
Strong institutions create structural immunity

Hurdle interpretation:

All countries could have coups, but the first coup requires extraordinary conditions
Once a country has one coup, subsequent coups become more likely (precedent exists)
Strong institutions make the first coup harder, not impossible

Which is right?

Both are plausible! The choice depends on your theory:

If you believe Switzerland is fundamentally different from Pakistan (qualitative difference), use ZIP
If you believe all countries face a threshold but some are much harder to cross (quantitative difference), use Hurdle
Empirically, compare models with Vuong test and information criteria
Substantively, consider which mechanism your theory emphasizes

6.4 Vuong Test for Non-Nested Models

6.4.1 Understanding Nested vs. Non-Nested Models

Before we can choose the right test for model comparison, we need to understand the concept of nesting.

Nested models: One model is a special case of another. You can get the simpler model by imposing constraints on the more complex model.

Examples of nested models:

Poisson vs. Negative Binomial: The Poisson is nested within the NB. If you set the overdispersion parameter $\alpha = 0$ in the NB model, you get the Poisson. The NB has all the Poisson parameters plus $\alpha$.
ZIP vs. ZINB: The ZIP is nested within the ZINB. Set $\alpha = 0$ in ZINB and you get ZIP.

When models are nested, we use the likelihood ratio test (LRT):

\[ \text{LR} = 2(\ell_{\text{complex}} - \ell_{\text{simple}}) \sim \chi^2_{\text{df}} \]

where df = difference in number of parameters. This test is valid because the simpler model is just the complex model with constraints.

Non-nested models: Neither model is a special case of the other. They have fundamentally different structures.

Examples of non-nested models:

ZIP vs. Hurdle: These are structurally different models:
- ZIP says: some zeros are structural (immune), others are sampling zeros
- Hurdle says: ALL zeros come from one process (failure to cross threshold)
You can’t get one from the other by setting parameters to special values. They’re different theoretical frameworks.
ZIP vs. Poisson: You might think ZIP nests Poisson (set zero-inflation to zero?), but it doesn’t work that way. When the zero-inflation probability is zero, ZIP has extra parameters that become meaningless. The models aren’t properly nested.
Hurdle vs. Poisson: Same issue. The Hurdle’s binary component doesn’t reduce to Poisson when adjusted.

6.4.2 Why Non-Nested Models Need Different Tests

The likelihood ratio test doesn’t work for non-nested models because:

Neither model is a constrained version of the other
We can’t say one model has “extra parameters” beyond the other
The LR statistic doesn’t follow a $\chi^2$ distribution

This is where the Vuong test comes in. Named after Quang Vuong (1989), it’s designed specifically for comparing non-nested models.

6.4.3 The Vuong Test: Intuition and Interpretation

The Vuong test asks: “Which model is closer to the true data-generating process?” It compares how well each model fits each observation, then tests whether one model systematically fits better.

The test statistic: For each observation $i$, compute the log-likelihood ratio:

\[ m_i = \log f_1(y_i | \hat{\boldsymbol{\theta}}_1) - \log f_2(y_i | \hat{\boldsymbol{\theta}}_2) \]

If Model 1 fits observation $i$ better, $m_i > 0$. If Model 2 fits better, $m_i < 0$. The Vuong statistic is:

\[ V = \frac{\sqrt{n} \times \bar{m}}{s_m} \]

where $\bar{m}$ is the average of the $m_i$ and $s_m$ is their standard deviation. This is essentially a t-test of whether the average log-likelihood difference is significantly different from zero.

Interpretation:

|V| < 1.96: Models fit equally well (at 5% level), can’t distinguish between them
V > 1.96: Model 1 fits significantly better than Model 2
V < -1.96: Model 2 fits significantly better than Model 1

Important caveats for political scientists:

The Vuong test is just one piece of evidence: Don’t choose models based solely on this test. Theory matters more.
Equal fit doesn’t mean both are good: If V is not significant, both models might fit poorly. Check absolute fit measures (rootograms, residual plots).
Sample size matters: With small samples (n < 100), the test has low power. With huge samples (n > 10,000), trivial differences become “significant.”
Not all Vuong tests are the same: Some software uses the original Vuong test, others use corrections (Wilson 2015, Merkle & You 2021). Check your package documentation.

6.4.4 Applying the Vuong Test

Let’s compare our ZIP and Poisson models:

# compare ZIP vs. standard Poisson (from earlier)
# note: we need to refit Poisson with same data as ZIP for fair comparison
vuong_result <- vuong(m_pois_zi, m_zip)

## Vuong Non-Nested Hypothesis Test-Statistic: 
## (test-statistic is asymptotically distributed N(0,1) under the
##  null that the models are indistinguishible)
## -------------------------------------------------------------
##               Vuong z-statistic             H_A    p-value
## Raw                   -13.53916 model2 > model1 < 2.22e-16
## AIC-corrected         -13.43954 model2 > model1 < 2.22e-16
## BIC-corrected         -13.16055 model2 > model1 < 2.22e-16

print(vuong_result)

## NULL

How to report Vuong test results:

Good: “The Vuong test favors the ZIP model over standard Poisson (V = 3.45, p < 0.001), suggesting that zero-inflation improves fit. However, substantively, we interpret this as evidence that some countries are structurally immune to protests given their institutional features.”

Bad: “The ZIP model is better (Vuong test significant).”

When to use Vuong vs. other tests:

Nested models (Poisson vs. NB, ZIP vs. ZINB): Use likelihood ratio test
Non-nested models (ZIP vs. Hurdle, ZIP vs. Poisson): Use Vuong test
Any models: Can always use AIC/BIC for comparison, though these are just information criteria, not formal tests

7 Model Selection and Practical Workflow

7.1 Decision Tree for Count Model Selection

A Step-by-Step Guide

7.1.1 Step 1: Is your outcome a true count?

YES: Proceed to Step 2
NO: Use appropriate model:
- Proportions (k of n) → Binomial/Beta regression
- Durations → Survival analysis
- Ordered categories → Ordered probit/logit
- Many categories (7+) → Consider OLS

7.1.2 Step 2: Check for dispersion

# quick check: compare mean and variance
mean(y); var(y)
# formal test after Poisson regression:
dispersiontest(poisson_model)

Var ≈ Mean: Use Poisson (proceed to Step 3)
Var > Mean: Use Negative Binomial (proceed to Step 3)
Var < Mean: Rare! Check for:
- Data errors
- Sample selection issues
- Consider underdispersed models (Conway-Maxwell-Poisson)

7.1.3 Step 3: Assess zero inflation

# compare observed vs predicted zeros
observed_zeros <- sum(y == 0)
predicted_zeros <- sum(dpois(0, fitted(model)))

Zeros match predictions: Stop here, use Poisson or NB
Excess zeros: Proceed to Step 4
Too few zeros: Check data quality, consider truncated models

7.1.4 Step 4: Choose between ZIP/ZINB and Hurdle

Ask yourself: “Why are there zeros?”

Structural immunity (some units can never experience outcome):

Switzerland can’t have coups (institutional impossibility)
Landlocked countries can’t have naval battles
→ Use Zero-Inflated models

Threshold/activation (all units could experience outcome, but need trigger):

All countries could protest, but need grievance + opportunity
All legislators could sponsor bills, but need motivation
→ Use Hurdle models

Not sure?

Fit both and compare with Vuong test
Consider theoretical implications of each
Check which makes better out-of-sample predictions

7.1.5 Step 5: Select covariates for each component

For Zero-Inflated models:

Zero-inflation part: What makes units immune?
Count part: What affects rate among at-risk units?

For Hurdle models:

Binary part: What triggers first event?
Count part: What affects frequency after first event?

Can use same or different covariates in each part

7.1.6 Quick Reference Table

Data Pattern	Model Choice	Key Assumption
Mean ≈ Variance, zeros as expected	Poisson	Events independent, constant rate
Overdispersion, zeros as expected	Negative Binomial	Unobserved heterogeneity
Overdispersion + excess zeros (immunity)	ZINB	Two types: immune vs at-risk
Overdispersion + excess zeros (threshold)	Hurdle NB	Different process for 0→1 vs 1→2
Underdispersion	Conway-Maxwell-Poisson	Competition/saturation effects

7.2 The Empirical Workflow

Start simple: Always begin with Poisson
- Establishes baseline
- Easiest to interpret
- May be adequate
Test dispersion: Use Cameron-Trivedi test
- If overdispersed → try negative binomial
- If underdispersed → rare, but consider quasi-Poisson
Check zero inflation: Compare observed vs predicted zeros
- Excess zeros → consider ZIP/ZINB or hurdle
- Theory should guide choice
Compare models: Use information criteria and tests
- Nested models: likelihood ratio test
- Non-nested: AIC/BIC, Vuong test
Validate: Check residuals and predictions
- Rootograms show fit across count distribution
- Residual plots reveal systematic problems

7.3 Common Pitfalls

Even experienced researchers make predictable errors with count models. Here are the most common mistakes and how to avoid them.

7.3.1 Pitfall 1: Ignoring Exposure

Comparing protest counts across countries without population adjustment is like comparing crime without considering city size. This is one of the most pervasive errors in applied count modeling.

The problem: You estimate a model predicting protest counts and find that India has more protests than Iceland. You conclude that India is more protest-prone. But India has 1.4 billion people while Iceland has 380,000. Even if both countries had identical per capita protest rates, India would have roughly 3,700 times more protests simply due to population size.

Why it matters: Without exposure adjustment, you’re confounding substance (political conditions that generate protests) with opportunity (population size that provides more potential protesters). Your coefficient on, say, GDP is picking up both “does GDP reduce protest propensity?” and “do richer countries have different population sizes?” You can’t separate these effects.

The fix: Always include exposure as an offset when appropriate:

Population exposure for events that scale with population (protests, crime, terrorism)
Temporal exposure for events that accumulate over time (conflicts during observation period)
Spatial exposure for events that scale with area (insurgent attacks across territory)

Example: Modeling protests without offset estimates how GDP affects raw protest counts (conflating population and propensity). Modeling with offset(log(population)) estimates how GDP affects protests per capita (the protest rate), which is usually what we care about theoretically.

Red flag: If your model finds that population size “doesn’t matter” when entered as a regular covariate, this often means you should be using it as an offset instead. The coefficient on population in count models is often close to 1, which is exactly what an offset assumes.

7.3.2 Pitfall 2: Zero-Inflating Everything

Not every excess zero needs a ZIP model. Sometimes you just need better covariates. Zero-inflation has become a methodological fashion that’s often misapplied.

The problem: You fit a Poisson model and observe that it underpredicts zeros. You immediately jump to ZIP or ZINB, get better fit statistics, and declare victory. But you haven’t asked whether the excess zeros indicate a genuine two-process mixture or just model misspecification.

Common causes of spurious excess zeros:

Omitted variables: You’re modeling terrorism but forgot to include state capacity. Strong states have zero terrorism not because they’re structurally immune, but because repression works. Adding state capacity as a covariate would reduce “excess” zeros without needing zero-inflation.
Wrong functional form: Democracy’s effect on protests might be non-linear (inverted U-shape). Forcing linearity creates excess zeros at the extremes (stable autocracies and consolidated democracies) that look like structural zeros but are just misspecified Poisson.
Ignored interactions: Oil wealth might suppress protests in autocracies but not democracies. Missing this interaction leaves systematic zeros in the “oil autocracy” group that look like zero-inflation.
Sample selection: Your data includes many country-years that shouldn’t be in the sample (e.g., countries during civil war when protest data isn’t collected). These are missing data zeros, not structural zeros.

The diagnostic: Before fitting ZIP/ZINB, try:

Adding theoretically relevant covariates you might have omitted
Testing polynomial terms or splines for key predictors
Including interactions that theory suggests
Examining who has zeros, are they a substantively distinct group?

If excess zeros persist after proper specification, then consider mixture models. But often, better specification eliminates the “excess” zeros.

Theoretical test: Ask yourself: “Do I have a theory about why some units are structurally immune?” If you can’t articulate such a theory, you probably don’t need zero-inflation. ZIP/ZINB are substantive models about dual processes, not just fixes for misspecification.

7.3.3 Pitfall 3: Misinterpreting ZIP Coefficients

ZIP coefficients are conditional on being at-risk. The overall marginal effect requires combining both parts. This is perhaps the most common interpretation error in published work.

The problem: You fit a ZIP model and report: “Democracy increases protests by 30% (exp(0.26) = 1.30).” But that 0.26 coefficient is from the count part of the model, which only applies to at-risk observations. You’ve ignored the zero-inflation part entirely, potentially getting the sign of the effect wrong.

Why it’s tricky: ZIP models produce two sets of coefficients:

Zero-inflation coefficients ($\boldsymbol{\gamma}$): Affect probability of being a structural zero
Count coefficients ($\boldsymbol{\beta}$): Affect the rate among at-risk units only

A variable can appear in both parts, potentially with opposite signs:

Democracy might decrease structural zero probability (makes protests possible) [negative $\gamma$]
But decrease count intensity among at-risk countries (institutionalizes dissent) [negative $\beta$]
Overall effect: Democracy increases protest likelihood but decreases intensity, you can’t tell this from either coefficient alone

The correct approach: Always report:

Zero-inflation interpretation: “What affects structural immunity?”
Count interpretation: “What affects intensity among at-risk units?”
Overall marginal effect: “What’s the net impact on expected counts?”

The overall marginal effect requires simulation or the delta method because it’s a non-linear combination of both sets of parameters. Don’t just report count coefficients and claim you’re done.

Example of getting it wrong: A study finds positive count coefficients for economic development and concludes “development increases conflict.” But the zero-inflation coefficients show development strongly predicts structural immunity. The overall effect might be negative, development reduces conflict by moving countries into the immune class, despite increasing conflict intensity among the few at-risk developed countries.

7.3.4 Pitfall 4: Forgetting About Theory

Statistical fit shouldn’t override theoretical sense. If your theory says all countries can experience protests (just at different rates), don’t use ZIP just because it has lower AIC. Model selection is a dialogue between theory and data, not a mechanical exercise.

The problem: You compare Poisson, NB, ZIP, and ZINB using AIC. ZINB wins. You use ZINB and move on. But you never asked whether zero-inflation makes theoretical sense for your application.

Why AIC isn’t enough: Information criteria measure statistical fit but ignore theoretical coherence. ZIP might fit better simply because it has more parameters (it’s more flexible), not because the two-process story is correct. With enough parameters, you can fit any data, but that doesn’t mean the model captures the true data-generating process.

Theoretical coherence questions:

For ZIP: Can you name specific units that are structurally immune? Or is zero-inflation just capturing unobserved heterogeneity that NB also handles?
For hurdle: Is crossing from 0 to 1 qualitatively different from 1 to 2? Or are you forcing a threshold that doesn’t exist?
For mixtures: Do you believe in distinct latent classes? Or is variation continuous?

Example of theory mattering: Modeling coups in OECD countries. ZIP might technically fit better (lots of zeros), but theoretically, consolidated democracies aren’t “structurally immune”, they just have extremely low rates. The zero-inflation interpretation (some countries CAN’T have coups) contradicts democratic theory (coups are always possible, just unlikely). Use NB with very low predicted rates instead.

The right approach:

Start with theory: What process generates the data?
Choose models consistent with that theory
Use diagnostics and fit statistics to select among theory-consistent models
If theory and fit conflict, investigate why, maybe theory is wrong, maybe data is limited

Red flag: If your best-fitting model contradicts your theoretical priors, don’t just accept the model. Either (a) you’ve discovered something theoretically important, or (b) the model is overfitting or misspecified. Investigate before concluding.

7.4 Complete Workflow Example

Let’s work through a full analysis of civil conflict onset:

# simulate civil war onset data
set.seed(6886)
n_countries <- 150
n_years <- 20
n <- n_countries * n_years

# country-year panel
country <- rep(1:n_countries, each = n_years)
year <- rep(1:n_years, n_countries)

# time-invariant country characteristics
ethnic_frac <- rep(runif(n_countries, 0, 1), each = n_years)
mountainous <- rep(runif(n_countries, 0, 100), each = n_years)

# time-varying characteristics
gdp_growth <- rnorm(n, 2, 3)
population <- exp(rnorm(n, 15, 1.5))  # log population
oil_price <- rep(rnorm(n_years, 50, 20), n_countries)  # global oil price (affects all countries)
democracy <- rbinom(n, 1, 0.3)

# generate conflict onsets
# theory: ethnic fractionalization and mountains increase risk
# economic growth reduces risk, democracy has complex effect
# high oil prices globally can increase conflict over resource wealth
beta0 <- -6
beta_ethnic <- 2
beta_mountain <- 0.02
beta_growth <- -0.15
beta_democracy <- 0.5
beta_oil <- 0.01

log_lambda <- beta0 + beta_ethnic*ethnic_frac + beta_mountain*mountainous +
              beta_growth*gdp_growth + beta_democracy*democracy +
              beta_oil*oil_price

lambda <- exp(log_lambda)

# add overdispersion and zero-inflation
# some countries are "immune" to civil war (structural zeros)
immune <- rep(rbinom(n_countries, 1, 0.2), each = n_years)  # 20% immune
conflicts <- ifelse(immune == 1, 0, rnbinom(n, size = 0.5, mu = lambda))

conflict_data <- data.frame(
  country = factor(country),
  year = year,
  conflicts = conflicts,
  ethnic_frac = ethnic_frac,
  mountainous = mountainous,
  gdp_growth = gdp_growth,
  population = population,
  democracy = democracy,
  oil_price = oil_price
)

# step 1: descriptive analysis
# descriptive statistics
desc_stats <- data.frame(
  Total_N = nrow(conflict_data),
  N_With_Conflict = sum(conflicts > 0),
  Pct_With_Conflict = round(mean(conflicts > 0) * 100, 1),
  Mean_All = round(mean(conflicts), 2),
  Mean_If_Positive = round(mean(conflicts[conflicts > 0]), 2),
  Variance = round(var(conflicts), 2),
  Var_Mean_Ratio = round(var(conflicts)/mean(conflicts), 2)
)
print(desc_stats)

##   Total_N N_With_Conflict Pct_With_Conflict Mean_All Mean_If_Positive Variance
## 1    3000             105               3.5     0.04             1.11     0.05
##   Var_Mean_Ratio
## 1            1.2

# fit progression of models
m1_pois <- glm(conflicts ~ ethnic_frac + mountainous + gdp_growth +
               democracy + oil_price + offset(log(population)),
               family = poisson, data = conflict_data)

test_result <- dispersion_test(m1_pois)

##              Delta     SE  t_stat p_value Var_Mean_Ratio
## lambda_hat -0.6158 2.8043 -0.2196  0.8262            1.2
##                                              Conclusion
## lambda_hat No significant dispersion (Poisson adequate)

# negative binomial (because overdispersion detected)
m2_nb <- glm.nb(conflicts ~ ethnic_frac + mountainous + gdp_growth +
                democracy + oil_price + offset(log(population)),
                data = conflict_data)
# check for excess zeros
zero_check <- data.frame(
  Observed_Zeros = sum(conflicts == 0),
  Expected_Zeros_NB = round(sum(dnbinom(0, size = m2_nb$theta, mu = fitted(m2_nb)))),
  Excess_Zeros = round(sum(conflicts == 0) - sum(dnbinom(0, size = m2_nb$theta, mu = fitted(m2_nb))))
)
print(zero_check)

##   Observed_Zeros Expected_Zeros_NB Excess_Zeros
## 1           2895              2900           -5

# ZIP, ZINB, and Hurdle models
m3_zip <- zeroinfl(conflicts ~ ethnic_frac + mountainous + gdp_growth +
                   democracy + oil_price + offset(log(population)) |
                   ethnic_frac + democracy,
                   dist = "poisson", data = conflict_data)

m4_zinb <- zeroinfl(conflicts ~ ethnic_frac + mountainous + gdp_growth +
                    democracy + oil_price + offset(log(population)) |
                    ethnic_frac + democracy,
                    dist = "negbin", data = conflict_data)

m5_hurdle <- hurdle(conflicts ~ ethnic_frac + mountainous + gdp_growth +
                    democracy + oil_price + offset(log(population)) |
                    ethnic_frac + democracy,
                    dist = "negbin", data = conflict_data)

# step 3: model comparison
# model comparison
model_list <- list(
  Poisson = m1_pois,
  NegBin = m2_nb,
  ZIP = m3_zip,
  ZINB = m4_zinb,
  Hurdle = m5_hurdle
)

comparison <- data.frame(
  Model = names(model_list),
  LogLik = round(sapply(model_list, logLik), 1),
  AIC = round(sapply(model_list, AIC), 1),
  BIC = round(sapply(model_list, BIC), 1)
)
print(comparison[order(comparison$AIC), ])

##           Model LogLik    AIC    BIC
## Hurdle   Hurdle -473.1  966.2 1026.2
## NegBin   NegBin -516.7 1047.5 1089.5
## ZINB       ZINB -516.5 1053.0 1113.1
## ZIP         ZIP -543.0 1104.1 1158.1
## Poisson Poisson -577.5 1167.0 1203.1

# best model coefficients (ZINB)
summary(m4_zinb)

## 
## Call:
## zeroinfl(formula = conflicts ~ ethnic_frac + mountainous + gdp_growth + 
##     democracy + oil_price + offset(log(population)) | ethnic_frac + democracy, 
##     data = conflict_data, dist = "negbin")
## 
## Pearson residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29455 -0.17515 -0.09863 -0.04719 43.42634 
## 
## Count model coefficients (negbin with log link):
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -21.093399   1.431473 -14.735  < 2e-16 ***
## ethnic_frac   2.907780   1.247601   2.331 0.019769 *  
## mountainous   0.027868   0.005364   5.195 2.04e-07 ***
## gdp_growth   -0.288402   0.055840  -5.165 2.41e-07 ***
## democracy     0.076977   0.620381   0.124 0.901252    
## oil_price    -0.001474   0.008128  -0.181 0.856089    
## Log(theta)   -2.309471   0.613308  -3.766 0.000166 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.01931    2.16051   0.009    0.993
## ethnic_frac -1.53951    2.15713  -0.714    0.475
## democracy   -0.62145    1.97536  -0.315    0.753
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Theta = 0.0993 
## Number of iterations in BFGS optimization: 35 
## Log-likelihood: -516.5 on 10 Df

count_results <- data.frame(
  Component = "Count",
  Variable = names(coef(m4_zinb, "count")),
  Coefficient = round(coef(m4_zinb, "count"), 3),
  IRR = round(exp(coef(m4_zinb, "count")), 3)
)

zero_results <- data.frame(
  Component = "Zero-Inflation",
  Variable = names(coef(m4_zinb, "zero")),
  Coefficient = round(coef(m4_zinb, "zero"), 3),
  Odds_Ratio = round(exp(coef(m4_zinb, "zero")), 3)
)

print(count_results)

##             Component    Variable Coefficient    IRR
## (Intercept)     Count (Intercept)     -21.093  0.000
## ethnic_frac     Count ethnic_frac       2.908 18.316
## mountainous     Count mountainous       0.028  1.028
## gdp_growth      Count  gdp_growth      -0.288  0.749
## democracy       Count   democracy       0.077  1.080
## oil_price       Count   oil_price      -0.001  0.999

print(zero_results)

##                  Component    Variable Coefficient Odds_Ratio
## (Intercept) Zero-Inflation (Intercept)       0.019      1.019
## ethnic_frac Zero-Inflation ethnic_frac      -1.540      0.214
## democracy   Zero-Inflation   democracy      -0.621      0.537

7.5 Visualization and Model Checking

A critical step is comparing predicted versus observed zeros, the key diagnostic for choosing between Poisson, NB, ZIP, and ZINB.

# compare predicted vs observed distribution
count_range <- 0:min(max(conflicts), 10)
observed_freq <- table(factor(conflicts, levels = count_range))

# get predicted probabilities from each model
pred_pois <- sapply(count_range, function(y) mean(dpois(y, fitted(m1_pois))))
pred_nb <- sapply(count_range, function(y) mean(dnbinom(y, size = m2_nb$theta, mu = fitted(m2_nb))))
pred_zinb <- colMeans(predict(m4_zinb, type = "prob")[, as.character(count_range)])

pred_df <- data.frame(
  Count = count_range,
  Observed = as.numeric(observed_freq),
  Poisson = round(pred_pois * nrow(conflict_data)),
  NegBin = round(pred_nb * nrow(conflict_data)),
  ZINB = round(pred_zinb * nrow(conflict_data))
)
print(pred_df)

##   Count Observed Poisson NegBin ZINB
## 0     0     2895    2904   2900 2900
## 1     1       95      82     55   54
## 2     2        8      10     16   17
## 3     3        2       3      8    8

The table shows observed counts versus predicted counts for each value. Models that severely underpredict zeros suggest need for zero-inflation or better covariates.

7.6 Reporting Results

When reporting count models:

Show multiple models: Typically main in the manuscript rest in the appendix
Report both statistics:
- Coefficients (log scale) with standard errors
- IRRs or percentage changes for interpretation
Substantive effects: Predicted counts for meaningful scenarios
Goodness of fit: AIC/BIC, dispersion tests, zero inflation tests
Robustness: Alternative specifications, sample restrictions

Example coefficient table across models:

key_models <- list("Poisson" = m1_pois, "Neg Binomial" = m2_nb, "ZINB" = m4_zinb)

ethnic_effects <- data.frame(
  Model = names(key_models),
  Coefficient = NA,
  SE = NA,
  IRR = NA
)

for(i in 1:length(key_models)) {
  m <- key_models[[i]]
  if(class(m)[1] %in% c("glm", "negbin")) {
    ethnic_effects$Coefficient[i] <- coef(m)["ethnic_frac"]
    ethnic_effects$SE[i] <- sqrt(vcov(m)["ethnic_frac", "ethnic_frac"])
  } else {
    ethnic_effects$Coefficient[i] <- coef(m, model = "count")["ethnic_frac"]
    ethnic_effects$SE[i] <- sqrt(vcov(m, model = "count")["ethnic_frac", "ethnic_frac"])
  }
  ethnic_effects$IRR[i] <- exp(ethnic_effects$Coefficient[i])
}

ethnic_effects[, 2:4] <- round(ethnic_effects[, 2:4], 3)
print(ethnic_effects)

##          Model Coefficient    SE    IRR
## 1      Poisson       2.513 0.349 12.337
## 2 Neg Binomial       3.439 0.461 31.164
## 3         ZINB       2.908 1.248 18.316

8 Wrapping up

Count models aren’t just Poisson regression. They’re a toolkit for discrete, non-negative outcomes with specific data-generating processes. The key insights:

Start simple, build up: Poisson → test dispersion → NB or ZIP/ZINB as needed
Theory guides specification: Zero-inflation vs hurdle depends on your substantive model
Interpretation matters: IRRs for main effects, careful with ZIP coefficients
Check everything: Dispersion, zeros, residuals, predictions
Substance over statistics: Better AIC doesn’t mean better science

We counts things that matter, wars, coups, protests, votes. Count models help us understand not just whether these events occur, but how often and why. Use them wisely.

8.1 Real Examples in Political Science Literature

Classic Applications of Count Models:

Poisson Models:
- King (1988, AJPS): “Statistical Models for Political Science Event Counts” - foundational piece
- Beck, Katz & Tucker (1998, AJPS): Binary time-series cross-section models
Negative Binomial Models:
- Gleditsch & Ward (2000, ISQ): Interstate dispute initiation
- Braithwaite (2010, JCR): Terrorism diffusion across countries
- Salehyan & Gleditsch (2006, IO): Refugee flows and conflict spread
Zero-Inflated Models:
- Bagozzi (2015, Political Analysis): Human rights violations with structural zeros
- Daxecker (2007, JCR): Electoral violence with immunity from violence
- Young (2019, JOP): Terrorist group formation and activity
Hurdle Models:
- Brandt et al. (2000, AJPS): Militarized interstate disputes
- Dietrich et al. (2019, APSR): Foreign aid allocation decisions
- Box-Steffensmeier et al. (2020, Political Analysis): Campaign contributions

Methodological Advances:

Cameron & Trivedi (2013): “Regression Analysis of Count Data” - comprehensive reference
Hilbe (2014): “Modeling Count Data” - practical guide
Zeileis et al. (2008, JSS): R implementation in pscl package

8.2 Key Formulas Reference

Poisson PMF: $P(Y=y|\lambda) = \frac{e^{-\lambda}\lambda^y}{y!}$

Negative Binomial Variance: $Var(Y) = \lambda + \alpha\lambda^2$

ZIP PMF: \[P(Y=0) = \psi + (1-\psi)e^{-\lambda}\] \[P(Y=y) = (1-\psi)\frac{e^{-\lambda}\lambda^y}{y!} \text{ for } y>0\]

IRR Interpretation: $\exp(\beta_k)$ = multiplicative change in expected count

Dispersion Test: Regress $u_i = \frac{(y_i-\hat{\lambda}_i)^2 - y_i}{\hat{\lambda}_i\sqrt{2}}$ on $\hat{\lambda}_i$

8.3 R Functions Quick Reference

# poisson regression
glm(y ~ x1 + x2, family = poisson, data = data)

# negative binomial
MASS::glm.nb(y ~ x1 + x2, data = data)

# zero-inflated Poisson
pscl::zeroinfl(y ~ x1 + x2 | z1 + z2, dist = "poisson", data = data)

# zero-inflated NB
pscl::zeroinfl(y ~ x1 + x2 | z1 + z2, dist = "negbin", data = data)

# hurdle model
pscl::hurdle(y ~ x1 + x2 | z1 + z2, dist = "poisson", data = data)

# dispersion test
AER::dispersiontest(model)

# model comparison
pscl::vuong(model1, model2)  # non-nested models
lrtest(model1, model2)        # nested models

# offsets for exposure
glm(y ~ x1 + x2 + offset(log(exposure)), family = poisson)

PLS 803: Count Models in Political Science

From Poisson to Zero-Inflated and Hurdle Models

Shahryar Minhas

2025-10-15