Event count models

Packages

Loading packages:

library(lme4)
library(dplyr)
library(stargazer)
library(ggplot2)
library(ggeffects)
library(nnet)

---

Data

Loading the data:

load('/Users/davidhamad/Documents/Cand.Scient.Pol/2. Semester/Statistical models beyond linear regression - applied statistics for political scientists/11-12) Count outcomes/df_yoshinaka.rda')

df <- df_yoshinaka %>% as.data.frame()

---

Exercise 1

Question 1: Estimate a poisson model and report the results in a table.

Answer 1: The Poisson model estimates how many legislative reports each MEP handles per month, using log(months) as an offset to account for time on the committee. The predictors are ideological distance to the EP group (the main hypothesis), distance to the national party (a competing explanation), and committee expertise (a control). All three are significant: both distance variables are negative, and experts get about 25% more reports.

# Estimate Poisson model 
mod.pois <- glm(nreports ~ 
                  nom_dist1_ep +
                  nom_dist1_nat +
                  ctee_expert +
                  offset(log(months)),
                family = "poisson",
                data = df)

# Report in a table
stargazer(mod.pois,
          title = "Poisson regression of legislative reports per MEP",
          dep.var.labels = "Number of reports",
          covariate.labels = c("Distance to EP group",
                               "Distance to national party",
                               "Committee expert"),
          digits = 2,
          type = "text")

Poisson regression of legislative reports per MEP
======================================================
                               Dependent variable:    
                           ---------------------------
                                Number of reports     
------------------------------------------------------
Distance to EP group                -0.85***          
                                     (0.28)           
                                                      
Distance to national party           -0.94**          
                                     (0.44)           
                                                      
Committee expert                     0.22***          
                                     (0.05)           
                                                      
Constant                            -3.08***          
                                     (0.04)           
                                                      
------------------------------------------------------
Observations                           898            
Log Likelihood                      -1,868.37         
Akaike Inf. Crit.                   3,744.74          
======================================================
Note:                      *p<0.1; **p<0.05; ***p<0.01

(exp(coef(mod.pois)["ctee_expert"]) - 1) * 100

ctee_expert 
   24.97171 

---

Question 2: Calculate and plot the standardized residuals of the model against the fitted/predicted values.

Answer 2: The plot shows the standardized residuals (z) against the fitted values from the Poisson model. Most residuals lie between approximately −2 and +2, but several observations have z-values well above +5, with the most extreme reaching about +13. The spread also increases with the fitted values, which suggests that the Poisson assumption of equidispersion may not hold.

# Calculate fitted values and standardized residuals
df.plot <- data.frame(
  fitted = predict(mod.pois, type = "response"),
  z      = residuals(mod.pois, type = "pearson")
)

# Plot
ggplot(df.plot, aes(x = fitted, y = z)) +
  geom_point(alpha = 0.4) +
  labs(x = "Fitted values",
       y = "Standardized residuals (z)",
       title = "Standardized residuals vs. fitted values") +
  theme_minimal()

---

Question 3: Draw a horizontal line at y = 0 and two additional horizontal lines at ±2 (or 1.96). The standardized residuals are supposed to be in this bound. How does it look? Can you explain the high z values?

Answer 3: Most standardized residuals fall within the ±2 bound, but a substantial number of observations lie well above +2, with several reaching values between +5 and +13. The high z-values indicate that the model substantially underestimates the number of reports for the most productive MEPs. This is a sign of overdispersion. The variance of the data exceeds what the Poisson model assumes. Likely sources include omitted variables (such as activity, chair status, or committee fixed effects) that distinguish highly productive MEPs from the average, or positive contagion (the tendency for MEPs who already have one report to receive more).

# Plot
ggplot(df.plot, aes(x = fitted, y = z)) +
  geom_point(alpha = 0.4) +
  geom_hline(yintercept = 0, color = "black") +
  geom_hline(yintercept = c(-2, 2), color = "red", linetype = "dashed") +
  labs(x = "Fitted values",
       y = "Standardized residuals (z)",
       title = "Standardized residuals vs. fitted values") +
  theme_minimal()

---

Question 4: Use the instructions in Gelman and Hill (2007), p. 114 to calculate a measure of overdispersion.

Answer 4: Following Gelman & Hill (2007), the estimated overdispersion is calculated as the sum of squared standardized residuals divided by the expected value n − k. The sum of squared residuals is 2,725.72, compared to the expected value of 894, giving an overdispersion ratio of 3.05. Since any value above 2 is considered large, this confirms substantial overdispersion in the data. The Poisson model’s standard errors are therefore too small, and the reported significance levels are overstated. This motivates fitting a model that allows for greater variance, such as a quasi-Poisson or negative binomial.

# Standardized residuals (Pearson)
z <- residuals(mod.pois, type = "pearson")

# Sum of squared standardized residuals
sum_z2 <- sum(z^2)

# Expected value (n - k)
n <- length(z)
k <- length(mod.pois$coefficients)

# Overdispersion ratio
ratio <- sum_z2 / (n - k)

# Print results (1)
sum_z2

[1] 2725.72

# Print results (2)
n - k

[1] 894

# Print results (2)
ratio

[1] 3.048904

---

Exercise 2

Question 1: Create an effect plot where you illustrate the effect of ideology on influence. Remember that the poisson backtransformation only implies taking the exponent!

Answer 1: The figure shows the predicted number of legislative reports as a function of ideological distance to the EP group median, holding the other predictors at their typical values. The y-axis is on the response scale, obtained by exponentiating the linear predictor. MEPs at the group median are predicted to receive about 2 reports, while those at the largest observed distance are predicted to receive less than 1. The relationship is non-linear due to the exponential link, with the effect being most pronounced at small distances and flattening out at larger ones. The shaded area represents the 95% confidence interval, which widens at extreme values where fewer observations exist.

ggpredict(mod.pois, terms = "nom_dist1_ep") |>
  plot() +
  labs(x = "Distance to EP group median",
       y = "Predicted number of reports",
       title = "Effect of ideological distance to EP group on legislative reports")