Outcome variable: number of legislative proposals

The outcome variable is the number of reports (i.e. legislative proposals, nreports) each MEP has garnered during the legislative period. We have, in other words, a data generating process in which events occur – in principle unrelated to each other – within a certain window/exposure time. This is the type of outcomes where we will resolve to the poisson model or one of its affiliates.

df %>%
  ggplot + 
  geom_histogram(
    aes(nreports)
    ) +
  xlab("Number of reports  (proposals) per MEPs") +
  ylab("Number of MEPs") +
  ggtitle("Distribution of legislative proposals among  MEPs")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

The distribution of legislative proposals is quite unequal. The distribution looks like a “ski slope”: I have many zero counts (MEPs who never proposed legislation) and some really high counts (MEPs that have handled quite a few proposals).

#Number of observations (MEPs)
nrow(df)

## [1] 906

#Number of no events (MEPs without any reports)
sum(df$nreports==0)

## [1] 367

#Number of extreme events
quantile(df$nreports, 0.95)

## 95% 
##   7

max(df$nreports)

## [1] 22

A fairly large number of MEPs (more than a third) never drafted a report at all, while the 5% most productive members drafted 7 or more reports. The most prolific MEP drafted no less than 22 proposals!

Overall, the variance in our outcome is higher than the mean.

#Mean
mean(df$nreports)

## [1] 1.970199

#Variance
var(df$nreports)

## [1] 8.340547

This is not uncommon for count data. The poisson distribution used in poisson models assume that the mean is the same as the spread/variance in the conditional outcome. If our regression model contains enough predictors to to explain this variance, that will not be a problem. However, it means I need predictors that explain both the long tail and the high zero-count.

Main predictors: political distance

Our research question concerns the political distance between each MEP and their political groups (i.e. the median “voter” in the group). Here, I explore the distribution of our main predictor nom_dist1_ep, the left-right orientation of each MEP and their absolute distance from their group. Out of curiosity, I plot it per political group.

df %>%
  ggplot + 
  geom_histogram(
    aes(nom_dist1_ep)
    ) +
  xlab("Left-rigth political distance between MEP and EP group median") +
  ylab("Number of MEPs") +
  ggtitle("Distribution of political distance among  MEPs") 

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 7 rows containing non-finite values (`stat_bin()`).

My interest in the distribution of the predictor is first and foremost to check what variation I have. I can see that there is quite some variation between MEPs, but I already get a sense that there are quite a few MEPs that are really close to their group, while some people are far off from the group median. It already tells me that – even if my hypothesis is correct – my predictor cannot explain all the variation in the dependent variable. That’s because I’d expect that the higher counts in reports would be among MEPs that are close to their leaders. Here, we see that the group leaders have many MEPs to choose from that are politically close to them; yet only very few MEPs are going to be “super-rapporteur” (i.e. long tails). Similarly, we have too many zero-counts in the outcome compared to political outliers in the predictor.

Control variables

We will include other control variables as well: The political distance to the national party, expertise, activity levels etc.

library(stargazer)
df_tab <- 
  df %>%
  #The select function from dplyr; not from MASS
  dplyr::select(nreports, nom_dist1_ep, nom_dist1_nat,
                ctee_expert, activity, attendance, months) %>%
  #Transform to a data frame for stargazer to "understand" (sometimes necessary)
  as.data.frame

stargazer(df_tab,
          title = "Descriptive statistics",
          #File type; could be "text"
          type = "html",
          #I want summary statistics from the data frame
          summary = T)

Why not a linear model? (what not to do)

Ward and Ahlquist (2018) warn us to not run linear regressions on count data. There is one exception to that: When all your predictors are categorical, then you’re just calculating the difference between group averages anyways; you could to a linear model.

In our case, however, we have a linear predictor and my base-line go-to is the poisson regression. The poisson distribution has a single paramter \(\lambda\) (lambda). It describes both its mean/expected value and the spread/variance. That’s both its biggest asset and a draw-back.

The asset is that the variance increases with the mean. That is, the distribution “knows” that your data will be heteroscedastic. Look:

mod1.ols <- lm(nreports ~ 
                 nom_dist1_ep +
                 nom_dist1_nat+
                 ctee_expert,
              df)

Standardized residuals:

#Predictions
df$preds_ols <- predict(mod1.ols, df)
#Residuals
df$resid_ols <- df$nreports - df$preds_ols
#Standardize to z-scores
df$resid_ols_z <- df$resid_ols/sd(df$resid_ols, na.rm = T)

Plot the residuals against the predicted values.

df %>%
  ggplot +
  geom_point(aes(y = resid_ols_z,
                 x = preds_ols)) +
  ggtitle("Residual plot")

## Warning: Removed 8 rows containing missing values (`geom_point()`).

We have a fan-shaped distribution. As my predictions increase, I will be increasingly wrong (i.e. my residuals increase). If I were to run an OLS, my standard errors would be too small, because the normal distribution assumes an equal spread.

Some people “fix” the problem by log-transforming the dependent variable before they run a linear model (drawing from a normal distribution). However, even though this fixes the problem that my distribution is bounded at zero, it does not fix the problem with heteroscedasticity.

mod2.ols <- lm(log(nreports+1) ~ 
                 nom_dist1_ep +
                 nom_dist1_nat+
                 ctee_expert,
              df)

Use the ready-made function to get a residual plot. Hit “Return” on your key-board to flip through different model statistics.

plot(mod2.ols)

The drawback of the assumption of an equal mean and spread, is that sometimes the spread is still too high, even for the poisson distribution. That’s when we have overdispersion.

Marginal effects

A one-unit change in x would change y by a factor of \(exp(\beta)\).

exp(mod2.pois$coefficients["nom_dist1_ep"])

## nom_dist1_ep 
##    0.4260668

Once I’ve taken the exponent, I see that the factor is less than one; i.e. the number of legislative proposals decrease as the political distance increases. I’m then doing my usual trick to change this into percentage change.

eff1 <- (1-exp(mod2.pois$coefficients["nom_dist1_ep"]))
eff1

## nom_dist1_ep 
##    0.5739332

A one-unit increase in the political distance decreases the number of legislative proposals handled by the MEP by 57%.

How likely is such a shift? That is, how much “elasticity” is there in my predictor?

To determine a likely shift in my x-variable – the political distance between the MEP and the group – I consider the interquartile range. That’s a shift across the 50% most likely values in the data.

#Consider
summary(df$nom_dist1_ep)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## 0.00000 0.01700 0.03900 0.07476 0.08200 1.20500       7

#Interquartile range (1st to 3rd quartiles)
iqr <- IQR(df$nom_dist1_ep, na.rm = T)
iqr

## [1] 0.06499999

Hmmm… a one unit change is extremely unlikely. In fact, a common distance is 0.065. Let’s reformulate, then.

eff2 <- 1-(exp(mod2.pois$coefficients["nom_dist1_ep"] * iqr))
eff2

## nom_dist1_ep 
##   0.05394572

Right, so: comparing an MEP that is fairly close to the group (1st quartile) to an MEP that is fairly far from the group median (3rd quartile), I find that the most distant MEP tends to write 5% fewer proposals on behalf of the group than the ideologically close MEP.

Another scenario, given the topic, could be to compare someone that shares the group’s views completely (x = 0) to a distant MEP (3rd quartile).

#Consider
summary(df$nom_dist1_ep)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## 0.00000 0.01700 0.03900 0.07476 0.08200 1.20500       7

#Interquartile range (1st to 3rd quartiles)
distance <- quantile(df$nom_dist1_ep, na.rm = T, probs = c(0.75))
distance

##   75% 
## 0.082

eff3 <- 1-(exp(mod2.pois$coefficients["nom_dist1_ep"] * distance))
eff3

## nom_dist1_ep 
##     0.067568

The relative effect is somewhat larger, but not massive.

First-differences

To calculate the first-difference, I’d have to set up some full-fledged scenarios, then fill in the equation for predictions, backtransform, then compare.

scenario <- data.frame(
  #Two scenarios: one close and one far away
  nom_dist1_ep = quantile(df$nom_dist1_ep, na.rm = T, probs = c(0, 0.75)),
  #Set the control to a typical value; here, the median
  nom_dist1_nat = median(df$nom_dist1_nat, na.rm = T),
  #A non-expert
  ctee_expert = 0,
  #NB: the offset variable also has to be in, but don't log transform!
  months = 60)

#The coefficients
a = mod2.pois$coefficients[1]
b1 = mod2.pois$coefficients[2]
b2 = mod2.pois$coefficients[3]
b3 = mod2.pois$coefficients[4]

#The equation; first scenario
pred1.hand <- 
  #intercept
  a * 1 +
  #first scenario, political distance == 0
  b1 * scenario[1,1]+
  #typical distance to national party
  b2 * scenario[1,2]+
  #committee expert
  b3 * scenario[1,3] +
  #The offset variable has a regression coefficient of 1
  1 * log(scenario[1,4])

#The equation; second scenario 
pred2.hand <- 
  #intercept
  a * 1 +
  #second scenario, political distance == 0
  b1 * scenario[2,1]+
  #typical distance to national party
  b2 * scenario[2,2]+
  #committee expert
  b3 * scenario[2,3]  +
  #offset
  1 * log(scenario[1,4])

#backtransform
exp(pred1.hand); 

## (Intercept) 
##     2.74206

exp(pred2.hand)

## (Intercept) 
##    2.556784

#first difference
exp(pred2.hand) - exp(pred1.hand)

## (Intercept) 
##  -0.1852755

predict(mod2.pois, scenario,  type = "response")

##       0%      75% 
## 2.742060 2.556784

preds1 <- predict(mod2.pois, scenario,  type = "response")
preds1

##       0%      75% 
## 2.742060 2.556784

Effectplot with the ggpredict

I can use the ggeffects package to plot the effect of political distance on influence. However, the package does not like it when we recode our variables in the equation. This is what I did when I set the offset. I therefore reestimate the model.

#A new offset variable
df$offset <- log(df$months)
#Re-estimate the model
mod2.pois <- glm(nreports ~ 
                   nom_dist1_ep +
                   nom_dist1_nat +
                   ctee_expert,
                 offset = offset,
                 family = "poisson",
              df)

I let the political distance slide across a reasonable spectrum.

library(ggeffects)
eff <- ggpredict(mod2.pois,
                 terms = list("nom_dist1_ep" = seq(0, 0.5, 0.01)))

Here is the “dirty” way of plotting.

eff %>% plot

R fits the y axis so that everything looks fantastic. I can build on this to create a better plot.

eff %>%
  ggplot +
  geom_ribbon(
    aes(
      ymax = conf.high,
      ymin = conf.low,
      x = x
    ),
    #Transparent ribbon of uncertainty
    alpha = 0.3
  ) +
  geom_line(
    aes(
      y = predicted,
      x = x
    )
  ) +
  ylim(0,3) +
  ylab("Number of legislative proposals") +
  xlab("Political distance from the leadership") +
  ggtitle("Effect of political distance on productivity (legislative proposals) among MEPs")

Histograms

I usually start the really intuitive way. How does the distribution of my predictions compare to the distribution of the observed counts?

#Predict outcomes
df$preds_pois2 <- predict(mod2.pois, 
                          newdata = df,
                          #Note that we should not use the argument name here; but we want exponentiated outcomes
                          "response")

#Plot two overlayed histograms
df %>%
  ggplot() +
  #Observed outcomes
  geom_histogram(aes(x = nreports),
                 #transparent color
                 alpha=0.5)+
  #Predicted outcomes
  geom_histogram(aes(x = preds_pois2),
                 #transparent red color
                 fill = "red", 
                 alpha =  0.3)+
  #ggtheme
  theme_minimal() +
  ggtitle("Distribution of observed vs. predicted outcomes") +
  ylab("Number of MEPs") +
  xlab("Number of reports (legislative proposals)")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 8 rows containing non-finite values (`stat_bin()`).

This looks dismal. We are predicting too few zero events and too few extreme counts. Clearly, the variation in the distribution of influence in the EP is larger than what our current model accounts for. It seems we are in the presence of overdispersion. It means that the standard errors are too small, since the poisson distribution has clear expectations about the variance.

Hanging rootogram

I can plot the same thing using hanging rootograms (Ward and Ahlquist 2018, p.). The countreg/topmodels package has a ready-made function for us. It is available for installation online.

#install.packages("distributions3"); hit "Yes" in R console
#install.packages("topmodels", repos="http://R-Forge.R-project.org")
library(topmodels)

#Rootogram
rootogram(mod1.pois,
          main = "Hanging rootogram of mod1.pois")

The plot gives us the same message: It uses the span of predicted values to create the bars, but compares within each predicted count the number of wrong predictions. Once again, we underpredict zeros as well as extreme values. That is, all bars that hang below the zero-line, are counts we underpredict. All bars that never reach the zero line are counts that are overpredicted.

Zero-inflated models

Here, I model a poisson count model. The package I use is the pscl package containing – among other things – functions for zero-inflated and hurdle models.

#install.packages("pscl")
library(pscl)

## Warning: pakke 'pscl' blev bygget under R version 4.3.1

I can also make a theory out of this. Let’s say that there is a group of MEPs that are never takers; they’d rather be caught dead than write a legislative proposal. That is, before the group leader makes the selection, there is an initial stage of self (de-)selection into legislative work.

One predictor for this is the activity level of the MEP. The variable attendance captures the share of committee meetings that the member has attended. The idea is that the less you attend a committee, the more you opt out of legislative work.

I specify the predictors for the two models in the function and separate them with a vertical bar “|”). Before the bar, I specify the predictors for the count part of the model. After the bar are the predictors for the zero/hurdle part (binomial logistic model). Here, I specify that Attendance is a predictor for whether an MEP is willing to engage in committee work, but not for how many positive counts will see.

mod1.zero.pois <- zeroinfl(nreports ~ 
                             nom_dist1_ep + 
                             nom_dist1_nat +
                             chair +
                             ctee_expert +
                             #Remove the fixed effects for convenience
                             # committee +
                             offset(log(months))|attendance,
                           #Count model == poisson
                           dist = "poisson",
                    df)

summary(mod1.zero.pois)

## 
## Call:
## zeroinfl(formula = nreports ~ nom_dist1_ep + nom_dist1_nat + chair + 
##     ctee_expert + offset(log(months)) | attendance, data = df, dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0975 -0.7359 -0.3825  0.3695  9.2827 
## 
## Count model coefficients (poisson with log link):
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.92656    0.03987 -73.408   <2e-16 ***
## nom_dist1_ep  -0.39240    0.36559  -1.073   0.2831    
## nom_dist1_nat -0.75077    0.43319  -1.733   0.0831 .  
## chair          1.32410    0.07129  18.573   <2e-16 ***
## ctee_expert    0.08714    0.05080   1.715   0.0863 .  
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.4276     0.2428   5.879 4.12e-09 ***
## attendance   -5.6132     0.5893  -9.525  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 14 
## Log-likelihood: -1548 on 7 Df

The effect comes out significant and with a massive substantial effect.

#Typical difference betweeen assiduous members and active members
change <- df$attendance %>%
  quantile(., c(0.2, .8), na.rm=T) %>% 
  diff %>%
  as.numeric()

change

## [1] 0.5028251

#Extract the relevant predictor from the count part of the model
coef <- mod1.zero.pois$coefficients$zero["attendance"]
effect <- exp(coef * change)-1
effect 

## attendance 
##  -0.940541

Comparing a low-attendee (i.e. low quantile) with an assiduous member (high quantile), I find a 94% relative decrease in an MEP’s likelihood of getting at least one legislative proposal from their leader.

mod2.zero.pois <- zeroinfl(nreports ~ 
                             nom_dist1_ep +
                             nom_dist1_nat +
                             ctee_expert +
                             chair +
                             attendance +
                             offset(log(months)),
                           #Count model == poisson
                           dist = "poisson",
                           df)

# summary(mod2.zero.pois)

library(pscl)
mod1.zero.nb <- zeroinfl(nreports ~ 
                      nom_dist1_ep +
                      nom_dist1_nat +
                      ctee_expert +
                        chair +
                        # committee +
                    offset(log(months)) | attendance,
                    #Count model == negative binomial
                    dist = "negbin",
                    df)

stargazer(mod1.zero.pois,mod2.zero.pois, mod1.zero.nb,
           omit = "committee",
          title = "The ZERO part of the zero-inflated models: poisson vs. the negative binomial model",
          dep.var.caption = "Zero events (no reports)",
          zero.component = T,
          type = "html")

stargazer(mod1.zero.pois,mod2.zero.pois, mod1.zero.nb,
           omit = "committee",
          title = "The COUNT PART of the zero-inflated models: poisson vs. the negative binomial model",
          dep.var.caption = "Number of reports",
          zero.component = F,
          type = "html")

Hurdle models

The hurdle models are similar. They assume that there is a certain threshhold an observation has to pass in order to acquire positive counts.

This time, we use the hurdle() function to estimate the model.

mod1.hurdle <- hurdle(nreports ~ 
                      nom_dist1_ep +
                      nom_dist1_nat +
                      ctee_expert +
                        chair +
                        # committee +
                    offset(log(months)) | attendance,
                    #Count model == negative binomial
                    # dist = "negbin",
                    df)

summary(mod1.hurdle)

## 
## Call:
## hurdle(formula = nreports ~ nom_dist1_ep + nom_dist1_nat + ctee_expert + 
##     chair + offset(log(months)) | attendance, data = df)
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3173 -0.7379 -0.4270  0.4003  8.7110 
## 
## Count model coefficients (truncated poisson with log link):
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.90541    0.04091 -71.011   <2e-16 ***
## nom_dist1_ep  -0.02383    0.36089  -0.066   0.9473    
## nom_dist1_nat -0.95820    0.52361  -1.830   0.0673 .  
## ctee_expert    0.06116    0.05263   1.162   0.2451    
## chair          1.29816    0.07142  18.177   <2e-16 ***
## Zero hurdle model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.5738     0.1771  -8.888   <2e-16 ***
## attendance    4.0378     0.3286  12.289   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 12 
## Log-likelihood: -1629 on 7 Df

As you can see from the regression table, these models often yield similar results.

stargazer(mod1.zero.pois,mod2.zero.pois, mod1.zero.nb, mod1.hurdle,
           omit = "committee",
          title = "The ZERO part of the zero-inflated models: poisson vs. the negative binomial model",
          dep.var.caption = "Zero events (no reports)",
          zero.component = T,
          type = "html")

stargazer(mod1.zero.pois, mod2.zero.pois, mod1.zero.nb, mod1.hurdle,
           omit = "committee",
          title = "The COUNT PART of the zero-inflated models: poisson vs. the negative binomial model",
          dep.var.caption = "Number of reports",
          zero.component = F,
          type = "html")

Assessing the zeroinflated and hurdle models

Assessing and presenting the results from the models of excess zero is somewhat harder because the predicted outcomes are a function of two processes/models. Which of those two should we focus on? Or should we do both?

The results of the binomial part (zero-part) reports the likelihood of positive/negative counts, and can be interpreted as such. However, the results of the count part of the model, should be interpreted as the effects of the covariates contitional on being in the group of MEPs eligible for reports (given the binomial model).

Predictions

The predict() function will provide predictions as per usual. However, since the binomial model moderates the results of the count part, we should think about which of the models we’d like to present.

The type= argument lets you decide which outcome you’re interested in. * “count” reports the results from the count part of the model condiitonal on passing the hurdle * “zero” reports the probability from the binomial model.

• “response” will yield the total expected outcome (with decimals). However these predictions are done at the MEP-level, and may therefore deviate quite substantially from the observed counts, due to the interaction of the two models.

To compare the model predictions, it might be more useful to compare the probabilities of the counts, i.e. the probability of receiving 0, 1, 2 etc legislative proposals from the group leaders.

• “prob” reports the probability of the count number of reports

This last option is interesting. Now, the function yields not a vector of predictions but a matrix. Each column corresponds to the count categories (the nreports), while the observations (MEPs) are reported in the rows. The cells report the probability for each MEP to be in a specific category of counts.

#Predict outcomes as a matrix
preds_zero<- predict(mod1.zero.pois, 
                        newdata = df,
                        #Specify the probability of a count
                        "prob")

dim(preds_zero)

## [1] 906  23

head(preds_zero)

##            0          1         2         3          4           5           6
## 1 0.23082468 0.11318799 0.1749544 0.1802843 0.13933247 0.086146324 0.044385377
## 2 0.31116564 0.09609989 0.1520713 0.1604280 0.12693297 0.080344966 0.042380070
## 3 0.08108013 0.11723650 0.1926832 0.2111220 0.17349404 0.114057945 0.062486367
## 4 0.36665427 0.29673189 0.2023309 0.0919748 0.03135717 0.008552533 0.001943889
## 5 0.30150344 0.26567402 0.2256174 0.1277335 0.05423735 0.018423918 0.005215360
## 6 0.08743681 0.12699462 0.2011780 0.2124635 0.16828654 0.106636155 0.056309067
##              7            8            9           10           11           12
## 1 0.0196018149 7.574618e-03 2.601792e-03 8.043166e-04 2.260418e-04 5.823203e-05
## 2 0.0191609874 7.580226e-03 2.665593e-03 8.436225e-04 2.427223e-04 6.401514e-05
## 3 0.0293425747 1.205644e-02 4.403395e-03 1.447434e-03 4.325307e-04 1.184805e-04
## 4 0.0003787053 6.455641e-05 9.781935e-06 1.333991e-06 1.653819e-07 1.879467e-08
## 5 0.0012654348 2.686602e-04 5.070075e-05 8.611284e-06 1.329624e-06 1.881920e-07
## 6 0.0254862179 1.009347e-02 3.553228e-03 1.125766e-03 3.242505e-04 8.561000e-05
##             13           14           15           16           17           18
## 1 1.384755e-05 3.057732e-06 6.301769e-07 1.217578e-07 2.214126e-08 3.802632e-09
## 2 1.558453e-05 3.523058e-06 7.433320e-07 1.470338e-07 2.737300e-08 4.812870e-09
## 3 2.995811e-05 7.033918e-06 1.541406e-06 3.166708e-07 6.123082e-08 1.118172e-08
## 4 1.971602e-09 1.920522e-10 1.746047e-11 1.488209e-12 1.193831e-13 9.044787e-15
## 5 2.458732e-08 2.982887e-09 3.377528e-10 3.585358e-11 3.582096e-12 3.380012e-13
## 6 2.086442e-05 4.721753e-06 9.973260e-07 1.974887e-07 3.680598e-08 6.478449e-09
##             19           20           21           22
## 1 6.187072e-10 9.563339e-11 1.407812e-11 1.978228e-12
## 2 8.016870e-10 1.268613e-10 1.911895e-11 2.750398e-12
## 3 1.934488e-09 3.179413e-10 4.976667e-11 7.435787e-12
## 4 6.491913e-16 4.426604e-17 2.874613e-18 1.781905e-19
## 5 3.021469e-14 2.565912e-15 2.075276e-16 1.602163e-17
## 6 1.080295e-09 1.711345e-10 2.581924e-11 3.718312e-12

This option will give us a better idea about the model’s performance in the aggregate because we can assess the distribution of the outcomes. We can either identify the most likely outcome as if each count were a category, or simply calculate the aggregated probabilities for each category/outcome.

Most likely outcome Here I sum over each column

colSums(preds_zero, na.rm = T) 

##            0            1            2            3            4            5 
## 328.49388994 128.11649395 132.15558796 111.25136593  78.76002226  47.94733459 
##            6            7            8            9           10           11 
##  25.70487224  12.60178980   6.06443721   3.23055103   2.13098970   1.71105210 
##           12           13           14           15           16           17 
##   1.49721003   1.31128497   1.10656961   0.88807610   0.67564952   0.48748810 
##           18           19           20           21           22 
##   0.33411593   0.21797666   0.13564867   0.08068526   0.04595930

which.max(preds_zero[1,]) 

## 0 
## 1

max(preds_zero[906,])  == preds_zero[906,]

##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE 
##    13    14    15    16    17    18    19    20    21    22 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

#Plot two overlayed histograms
df %>%
  ggplot() +
  #Observed outcomes
  geom_histogram(aes(x = nreports),
                 #transparent color
                 alpha=0.5)+
  #Predicted outcomes
  geom_histogram(aes(x = preds_pois2),
                 #transparent red color
                 fill = "red", 
                 alpha =  0.3)+
  #ggtheme
  theme_minimal() +
  ggtitle("Distribution of observed vs. predicted outcomes") +
  ylab("Number of MEPs") +
  xlab("Number of reports (legislative proposals)")

Here, we reported the most likely count category for each MEP. As is clear from the overlay, my model doesn’t predict the categories very precisely. That is, the probability of an MEP to fall into a single count category of reports is not very precisely estimated. However, if we treat the outcome as metric rather than categorical, the perception might change.

Distribution of probabilities We can simply compare the proportion of observed data in each category of nreports with the predicted probabilities. That way, we can take into account that an MEP has a fair chance of getting – say – 3-5 reports – but a low probability of getting 0 or 10.

To do so, I create two data sets – one with predicted values and one with observed values – then merge them. Last, I make a histogram.

#Data frame with predicted probabilities
predicted_prob <- 
  #the matrix from the predict function
  preds_zero %>%
  #make it into a data frame
  as.data.frame %>%
  #Take the mean of all the columns
  reframe(across(everything(), mean, na.rm = T)) %>%
  #Flip columns and rows to get a data set with two variables
  tidyr::pivot_longer(everything(),
                      names_to = "nreports",
                      values_to = "pred_prob") %>%
  #Recode to a numeric variable
  mutate(nreports = as.numeric(nreports))

#check
head(predicted_prob)

#Data frame with observed proportions
observed_prob <- 
  df %>%
  #Count the number of cases in each category
  group_by(nreports) %>%
  summarise(count = n()) %>%
  #Ungroup to consider all the data
  ungroup %>%
  #Calculate the number of each category
  mutate(N = sum(count),
         #divide by all observations to get a proportion
         prob = count/N,
         #Recode count category to numeric
         nreports = as.numeric(nreports))


#Merge the two data frames
df_probs <- full_join(observed_prob, 
                      predicted_prob,
                      by = "nreports")

#Plot
df_probs %>%
  # You can filter away the zeros to zoom in on the counts
  # filter(nreports > 0 ) %>%
  ggplot +
  geom_histogram(aes(y = prob,
                     x = nreports),
                 stat = "identity") +
  #Red == predicted values
  geom_histogram(aes(y = pred_prob,
                     x = nreports),
                 stat = "identity",
                 fill = "red", alpha = 0.5) +
  #Cut the x-axis to the values you want to see
  coord_cartesian(xlim = c(0,25)) +
  ggtitle("Predicted vs. observed counts of legislative proposals",
          subtitle = "Results from a zeroinflated model") +
  ylab("Probability") +
  xlab("Number of legislative proposals (`reports`)") +
  theme_minimal()

rootogram(mod1.zero.pois)