# Load relevant libraries
library(gbm) # gradient boosting machines: the models that predict AP and TDE
library(fixest) # fast fixed-effects estimation; supplies fenegbin() below
# plyr is NOT attached. Its ddply() is called as plyr::ddply() below.
# Attaching plyr would mask dplyr's summarise(), mutate(), rename() and
# arrange() -- plyr has functions of the same names that ignore dplyr
# grouping -- and that masking would persist for any other chapter
# rendered afterwards in the same R session.
library(tidyverse) # a bundle of data-handling packages (dplyr, ggplot2, ...)
library(mice) # multiple imputation; used here to pool results across datasets
library(miscTools) # assorted utilities; supplies colMedians()
library(kableExtra) # formats knitr tables for HTML output

### Function used to create expected values and SEs
# Adapted from https://www.r-bloggers.com/2013/01/the-cluster-bootstrap/

# ---------------------------------------------------------------------------
# ANATOMY OF A FUNCTION
# The four names in parentheses are the function's arguments -- the dials the
# user can turn. Those written with "= something" have a default value and may
# be omitted when the function is called.
#
#   myvalueCLPR  the level of personalism (AP) under the OLD closed-list system
#   myvalueFLPR  the level of personalism (AP) under the NEW free-list system
#   reps         how many bootstrap resamples to draw (500 by default)
#   seed         a fixed starting point for the random number generator
#
# Everything between the outermost { and } is the body: the instructions R
# carries out each time the function is called.
# ---------------------------------------------------------------------------
boot.predict.pork <- function(myvalueCLPR = 0.5, myvalueFLPR, reps = 500, seed = 123354) {
  # require() loads a package from inside a function, so the function still
  # works if it is copied into a fresh session where these are not yet loaded.
  require(fixest)
  require(miscTools)

  # Set Seed
  # Bootstrapping draws random samples. Fixing the seed means the "random"
  # draws are the same every time this code is run, so the numbers reported in
  # the book can be reproduced exactly.
  set.seed(seed)

  # Run Model
  # An empty list that will hold the five fitted models -- one per imputed
  # dataset. Creating the container first, then filling it, is standard R
  # practice.
  temp_list <- list()
  temp_list[[1]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[1]])
  temp_list[[2]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[2]])
  temp_list[[3]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[3]])
  temp_list[[4]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[4]])
  temp_list[[5]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[5]])

  # The five lines above fit the SAME model to the five imputed datasets.
  # Reading one of them from left to right:
  #   fenegbin()   fits a negative binomial regression. The outcome is a count
  #                (number of bills), which is why an ordinary linear model
  #                would be inappropriate: counts cannot go below zero and
  #                their variance grows with their mean.
  #   lpgbills ~   the ~ separates the outcome (left) from the predictors
  #                (right); read it as "is explained by".
  #   pers.hat     predicted Average Personalism -- the variable of interest
  #   totalEff.hat predicted Total Duvergerian Effect (interparty competition)
  #   govdeps ... sex  legislator-level controls
  #   se = 'hetero'    report heteroskedasticity-robust standard errors
  #   data =           which of the five datasets to use

  # smalldata
  # Keep just the predictor columns from the first dataset. The square brackets
  # index the data frame as [rows, columns]; leaving the row slot empty before
  # the comma means "all rows". c(...) lists the columns to retain.
  smallData <- hond_list[[1]][, c("pers.hat", "totalEff.hat", "govdeps", "seniority", "boardm", "permcomm", "smallparties", "sex")]

  ##### CLPR
  # Point Estimate
  # Build a single synthetic "typical legislator": the median value of every
  # predictor. colMedians() returns those medians as a row of numbers, t()
  # transposes it into one row, and as.data.frame() makes it something the
  # predict() function will accept.
  newData <- as.data.frame(t(miscTools::colMedians(smallData))) # copy data
  # Now override just one field: set personalism to its PRE-reform level and
  # leave every other characteristic at its median. This is the standard way of
  # isolating the effect of a single variable.
  # (The trailing comment is inherited from the blog post this was adapted
  # from, where the variable of interest was education rather than AP.)
  newData$pers.hat <- myvalueCLPR # Replace education value
  # Ask each of the five fitted models how many bills it expects for that
  # legislator. type = 'response' returns the answer on the scale of the
  # outcome -- an expected count of bills -- rather than on the model's
  # internal logarithmic scale.
  point.estimate1 <- predict(temp_list[[1]], newdata = newData, type = "response")
  point.estimate2 <- predict(temp_list[[2]], newdata = newData, type = "response")
  point.estimate3 <- predict(temp_list[[3]], newdata = newData, type = "response")
  point.estimate4 <- predict(temp_list[[4]], newdata = newData, type = "response")
  point.estimate5 <- predict(temp_list[[5]], newdata = newData, type = "response")
  # Avg
  # Rubin's first rule: the pooled point estimate is simply the average of the
  # five estimates, one from each imputed dataset.
  point.estimateCLPR <- (point.estimate1 + point.estimate2 + point.estimate3 + point.estimate4 + point.estimate5) / 5

  ##### FLrPR
  # Point Estimate
  # Repeat the identical exercise, changing only personalism to its POST-reform
  # level. Because newData is otherwise untouched, the difference between the
  # two answers is attributable to personalism alone.
  newData$pers.hat <- myvalueFLPR # Replace education value
  point.estimate1 <- predict(temp_list[[1]], newdata = newData, type = "response")
  point.estimate2 <- predict(temp_list[[2]], newdata = newData, type = "response")
  point.estimate3 <- predict(temp_list[[3]], newdata = newData, type = "response")
  point.estimate4 <- predict(temp_list[[4]], newdata = newData, type = "response")
  point.estimate5 <- predict(temp_list[[5]], newdata = newData, type = "response")
  # Avg
  point.estimateFLPR <- (point.estimate1 + point.estimate2 + point.estimate3 + point.estimate4 + point.estimate5) / 5

  ##### Difference
  # The quantity the chapter reports: how many additional local-goods bills the
  # typical legislator is expected to introduce once personalism rises to its
  # post-reform level.
  point.estimateDiff <- point.estimateFLPR - point.estimateCLPR

  # Calculate SE
  # ------------------------------------------------------------------------
  # THE BOOTSTRAP
  # A point estimate on its own says nothing about precision. The bootstrap
  # answers "how much would this number move if we had drawn a different
  # sample?" by repeatedly drawing new samples FROM THE DATA WE HAVE, with
  # replacement, and re-running the whole calculation on each one. The spread
  # of the resulting values estimates the standard error.
  #
  # Three storage matrices, each with one row per resample and one column per
  # imputed dataset. matrix(NA, ...) pre-fills them with missing values, which
  # get overwritten as the loop proceeds.
  # ------------------------------------------------------------------------
  sterrsCLPR <- matrix(NA, nrow = reps, ncol = 5)
  sterrsFLPR <- matrix(NA, nrow = reps, ncol = 5)
  sterrsDiff <- matrix(NA, nrow = reps, ncol = 5)
  # for(i in 1:reps) repeats everything inside the braces `reps` times, with i
  # taking the values 1, 2, 3, ... in turn.
  for (i in 1:reps) {
    # Create Data
    # Draw a resample: pick row numbers at random, as many as the dataset has
    # rows, with replacement (so some rows appear twice and others not at all).
    # The SAME row numbers are then applied to all five imputed datasets, which
    # keeps them aligned with one another.
    index <- sample(1:nrow(hond_list[[1]]), nrow(hond_list[[1]]), replace = TRUE) # sample clusters
    # Run Model
    # Re-fit the same five models on the resampled data. `hond_list[[1]][index, ]`
    # means "take these rows, all columns".
    temp_list_boot <- list()
    temp_list_boot[[1]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[1]][index, ])
    temp_list_boot[[2]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[2]][index, ])
    temp_list_boot[[3]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[3]][index, ])
    temp_list_boot[[4]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[4]][index, ])
    temp_list_boot[[5]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[5]][index, ])

    ### CLPR
    # Calculate predict value
    # Rebuild the "typical legislator" from the RESAMPLED data, so that
    # uncertainty about the medians is carried through as well.
    newDataBootData <- as.data.frame(t(miscTools::colMedians(smallData[index, ]))) # copy data
    newDataBootData$pers.hat <- myvalueCLPR # Replace education value
    # Store this resample's five predictions in row i of the matrix. `[i, 1]`
    # means row i, column 1.
    sterrsCLPR[i, 1] <- predict(temp_list_boot[[1]], newdata = newDataBootData, type = "response")
    sterrsCLPR[i, 2] <- predict(temp_list_boot[[2]], newdata = newDataBootData, type = "response")
    sterrsCLPR[i, 3] <- predict(temp_list_boot[[3]], newdata = newDataBootData, type = "response")
    sterrsCLPR[i, 4] <- predict(temp_list_boot[[4]], newdata = newDataBootData, type = "response")
    sterrsCLPR[i, 5] <- predict(temp_list_boot[[5]], newdata = newDataBootData, type = "response")

    ### FLPR
    # Calculate predict value
    # Same synthetic legislator, personalism switched to its post-reform level.
    newDataBootData$pers.hat <- myvalueFLPR # Replace education value
    sterrsFLPR[i, 1] <- predict(temp_list_boot[[1]], newdata = newDataBootData, type = "response")
    sterrsFLPR[i, 2] <- predict(temp_list_boot[[2]], newdata = newDataBootData, type = "response")
    sterrsFLPR[i, 3] <- predict(temp_list_boot[[3]], newdata = newDataBootData, type = "response")
    sterrsFLPR[i, 4] <- predict(temp_list_boot[[4]], newdata = newDataBootData, type = "response")
    sterrsFLPR[i, 5] <- predict(temp_list_boot[[5]], newdata = newDataBootData, type = "response")

    ### Difference
    # Column-by-column difference between the post- and pre-reform predictions.
    # Omitting the row index (`[, 1]`) operates on the whole column at once --
    # R's vectorised arithmetic, which removes the need for an inner loop.
    sterrsDiff[, 1] <- sterrsFLPR[, 1] - sterrsCLPR[, 1]
    sterrsDiff[, 2] <- sterrsFLPR[, 2] - sterrsCLPR[, 2]
    sterrsDiff[, 3] <- sterrsFLPR[, 3] - sterrsCLPR[, 3]
    sterrsDiff[, 4] <- sterrsFLPR[, 4] - sterrsCLPR[, 4]
    sterrsDiff[, 5] <- sterrsFLPR[, 5] - sterrsCLPR[, 5]
  }
  # prepare to return results
  # ------------------------------------------------------------------------
  # RUBIN'S SECOND RULE
  # With imputed data, total uncertainty has two sources, and the three lines
  # below add them together:
  #   WITHIN-imputation variance -- mean(apply(..., 2, var)) computes the
  #     variance down each of the five columns and averages them. This is the
  #     ordinary sampling uncertainty we would face even with complete data.
  #   BETWEEN-imputation variance -- the sum-of-squares term measures how far
  #     the five column means sit from the overall mean. It captures the extra
  #     uncertainty created by not knowing the missing values. The multiplier
  #     (5+1)/(5*(5-1)) is Rubin's correction for having only five imputations.
  # apply(X, 2, var) means "apply var() to each column of X" (2 = columns,
  # 1 would mean rows). na.rm = TRUE tells R to ignore missing values.
  # ------------------------------------------------------------------------
  sd_pool_CLPR <- mean(apply(sterrsCLPR, 2, var, na.rm = TRUE), na.rm = TRUE) + (5 + 1) / (5 * (5 - 1)) * (sum((colMeans(sterrsCLPR, na.rm = TRUE) - mean(sterrsCLPR, na.rm = TRUE))^2))
  sd_pool_FLPR <- mean(apply(sterrsFLPR, 2, var, na.rm = TRUE), na.rm = TRUE) + (5 + 1) / (5 * (5 - 1)) * (sum((colMeans(sterrsFLPR, na.rm = TRUE) - mean(sterrsFLPR, na.rm = TRUE))^2))
  sd_pool_Diff <- mean(apply(sterrsDiff, 2, var, na.rm = TRUE), na.rm = TRUE) + (5 + 1) / (5 * (5 - 1)) * (sum((colMeans(sterrsDiff, na.rm = TRUE) - mean(sterrsDiff, na.rm = TRUE))^2))

  # Bundle the six results into a single named vector. c() combines values;
  # names() labels them so the output is self-documenting when printed.
  val <- c(point.estimateCLPR, sd_pool_CLPR, point.estimateFLPR, sd_pool_FLPR, point.estimateDiff, sd_pool_Diff)
  names(val) <- c("Point Estimate - CLPR", "Std. Error - CLPR", "Point Estimate - FLPR", "Std. Error - FLPR", "Point Estimate - Diff", "Std. Error - Diff")
  # return() hands this vector back to whoever called the function.
  return(val)
}

# Function to extract estimates for table 7.4
extract_Estimates <- function(model_name) {
  # -------------------------------------------------------------------------
  # ifelse(test, value_if_true, value_if_false) is a vectorised choice: it is
  # applied to every coefficient at once. Here three ifelse() calls are nested
  # inside one another to walk down the conventional significance thresholds:
  #
  #   p < 0.01  ->  three stars
  #   p < 0.05  ->  two stars      (only reached if the first test failed)
  #   p < 0.10  ->  one star
  #   otherwise ->  no stars
  #
  # paste0() glues pieces of text together with no separator, so
  # paste0(0.123, "***", " (", 0.045, ")") produces "0.123*** (0.045)".
  # round(x, 3) keeps three decimal places.
  # -------------------------------------------------------------------------
  out <- ifelse(model_name$p.value < 0.01, paste0(round(model_name$estimate, 3), "***", " (", round(model_name$std.error, 3), ")"),
    ifelse(model_name$p.value < 0.05, paste0(round(model_name$estimate, 3), "**", " (", round(model_name$std.error, 3), ")"),
      ifelse(model_name$p.value < 0.1, paste0(round(model_name$estimate, 3), "*", " (", round(model_name$std.error, 3), ")"),
        paste0(round(model_name$estimate, 3), " (", round(model_name$std.error, 3), ")")
      )
    )
  )
  # Label each formatted string with the name of the variable it belongs to.
  names(out) <- model_name$term
  # Convert to a data frame, then copy the row names into an ordinary column
  # called "variable". That column is what merge() will match on later, when
  # the four models are combined side by side into one table.
  out <- data.frame(out)
  out$variable <- row.names(out)
  return(out)
}

# Set the path to Chapter 14 materials
# A line starting with # is inert: R skips it entirely. The line below is a
# switched-off duplicate of the active one, left over from a machine where the
# materials sat elsewhere. It has no effect and can be deleted.

# so the shorter file names below resolve correctly. This is the one line you
# must edit to run the chapter on your own machine.

# Load Honduras data
# load() reads a .RData file and restores whatever objects were saved in it,
# under their original names. Unlike readRDS(), you do not assign the result to
# anything -- the objects simply appear. This one supplies `hond`, the
# legislator-level dataset from Munoz-Portillo (2021).
load(file = "data/ch14/RData/hondData.RData")

# Load GBM objects
# These two supply `totalAP.objects` and `totalEffENP.objects`: the fitted
# gradient boosting models for intraparty (AP) and interparty (TDE) incentives.
load(file = "data/shared/AP_district_objects_t285_d13.RData")
load(file = "data/shared/TDE_district_objects_t485_d9.RData")

### Load GBM models used to predict AP and TDE

# Model for intraparty competition (AP)
# The $ operator pulls one named component out of a list. Each of these
# "objects" bundles were saved in Chapter 6 and contains, among other things,
# the best-performing fitted model under the name optimalGBM.
optimalGBMIntra <- totalAP.objects$optimalGBM

# Model for interparty dispersion (TDE)
optimalGBMInter <- totalEffENP.objects$optimalGBM

# Remove object to save memory
# The full bundles are large and no longer needed now that the models have been
# extracted. rm() deletes them from memory.
rm(totalAP.objects)
rm(totalEffENP.objects)

# Define all district magnitudes used in Honduras
# M = district magnitude, the number of seats a district elects. Honduras has
# districts of these eleven sizes, from one seat up to twenty-three.
M <- c(3, 4, 9, 23, 20, 7, 5, 6, 8, 1, 2)

# Create Honduras Dataset
# Build a description of the electoral system, one row per magnitude per era.
#   data.frame() creates a table; because M has eleven values while the other
#   arguments have one, R recycles the single values down all eleven rows.
#   rbind() ("row bind") stacks the two tables on top of each other, giving
#   22 rows in total.
#   tballot = 1 marks the POST-reform free-list system (open ballot);
#   tballot = 0 marks the PRE-reform closed-list system.
# The remaining columns are the other component rules the GBMs expect as
# inputs: the seat allocation formula, how votes are counted, the level at
# which votes are pooled, and any legal threshold.
hon <- rbind(
  data.frame(
    tballot = 1, ballot_type = "open", M = M, formula = "hare", new.nvotes =
      "TotalSeats", pool_level = "party", threshold = 0
  ),
  data.frame(
    tballot = 0, ballot_type = "closed", M = M, formula = "hare", new.nvotes = "One", pool_level =
      "party", threshold = 0
  )
)


# Adjust values for the single-member district (M == 1), where plurality applies
# A one-seat district cannot have an open list or a proportional formula, so
# those rows are corrected regardless of era. Read `hon$ballot_type[hon$M == 1]`
# as: "the ballot_type entries of hon, for those rows where M equals 1".
# Note the double equals sign: == asks a question, = assigns a value.
hon$ballot_type[hon$M == 1] <- "closed"
hon$formula[hon$M == 1] <- "plurality"
hon$new.nvotes[hon$M == 1] <- "One"

# Ensure categorical variables match the factor levels used in the GBM models
# A "factor" is R's type for categorical data: the labels are stored alongside
# a fixed list of permitted values, called levels. A fitted model remembers the
# levels it was trained on, and will refuse to predict -- or worse, will
# silently mismatch categories -- if new data present them in a different order.
# These four lines therefore copy the training levels straight out of the model
# object (`$var.levels`) and impose them on the new data.
hon$new.nvotes <- factor(hon$new.nvotes, levels = optimalGBMIntra[[1]]$var.levels[[3]])
hon$pool_level <- factor(hon$pool_level, levels = optimalGBMIntra[[1]]$var.levels[[2]])
hon$ballot_type <- factor(hon$ballot_type, levels = optimalGBMIntra[[1]]$var.levels[[1]])
hon$formula <- factor(hon$formula, levels = optimalGBMIntra[[1]]$var.levels[[6]])


# Reorder the data
# Sort rows by ballot type, then by magnitude within ballot type. order()
# returns the row numbers in sorted sequence; using them as a row index
# rearranges the data frame. The trailing comma means "keep all columns".
# Sorting here matters because later code pairs the rows up positionally.
hon <- hon[order(hon$ballot_type, hon$M), ]

### We need to run predictions five times
# Predict Interparty and Intraparty
# Five GBMs were trained (one per imputed dataset), so five sets of predictions
# are produced and stored in a list.
honData <- list()
for (i in 1:5) {
  # Interparty
  # predict() applies model i to the 22-row description of Honduras.
  # n.trees tells the boosting model how many trees to use -- here, all of the
  # trees selected as optimal during training.
  hon$totalEff.hat <- predict(optimalGBMInter[[i]], hon, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  hon$pers.hat <- predict(optimalGBMIntra[[i]], hon, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  # Store this iteration's copy of `hon`, predictions included, as element i.
  honData[[i]] <- hon
}

# Get TDE and AP, then average the values across datasets
# sapply(1:5, function(i) ...) runs the little unnamed function once for each
# i, extracting four columns from each of the five prediction sets;
# simplify = FALSE keeps the result as a list. do.call('rbind', ...) then
# stacks those five tables into one tall table of 5 x 22 = 110 rows.
all_TDE_AP <- do.call("rbind", sapply(1:5, function(i) honData[[i]][, c("tballot", "M", "totalEff.hat", "pers.hat")], simplify = FALSE))
# ddply() splits a data frame into groups, applies a function to each, and
# recombines the results. Here the groups are every combination of ballot era
# and district magnitude, and the function averages the five predictions --
# giving one pooled AP and one pooled TDE per system-magnitude pair.
avg_TDE_AP <- plyr::ddply(
  .data = all_TDE_AP, .variables = c("tballot", "M"), .fun = plyr::summarise,
  totalEff.hat = mean(totalEff.hat),
  pers.hat = mean(pers.hat)
)

# District-Year
# From here on we work with `hond`, the legislator-level data. paste0() glues
# two columns together to create a combined identifier, e.g. "Cortes-1997".
hond$district_year <- paste0(hond$district, "-", hond$year)

# Reform dummy
# hond$year > 2003 gives TRUE/FALSE; as.numeric() converts that to 1/0, the
# form regression models expect.
hond$reform <- as.numeric(hond$year > 2003)

# District-Reform
hond$district_reform <- paste0(hond$district, "-", hond$reform)


# Add AP and TDE to replication data
# Attach the predicted AP and TDE to each legislator-year observation.
hond_list <- list()
for (i in 1:5) {
  # Drop the unused factor levels by round-tripping through character. Levels
  # that appear in the model's training data but not in Honduras would
  # otherwise linger and confuse the merge.
  honData[[i]]$ballot_type <- as.factor(as.character(honData[[i]]$ballot_type))
  # Merge
  # merge() joins two tables wherever the matching columns agree. Here the keys
  # are ballot era (tballot) and district magnitude (M): every legislator gets
  # the AP and TDE values belonging to the system they were elected under.
  # The result is the list of five analysis datasets used by every model below.
  hond_list[[i]] <- merge(hond, honData[[i]][, c("tballot", "M", "totalEff.hat", "pers.hat", "ballot_type")], by.x = c("tballot", "M"))
}

## Estimate four nested negative binomial models across 5 imputed datasets
# "Nested" means each model contains all the terms of the one before it, plus
# something extra. Four empty lists, one per specification; each will hold five
# fitted models.
# Create empty lists to store models
m1_list <- list() # Model 1: AP only
m2_list <- list() # Model 2: AP + TDE
m3_list <- list() # Model 3: AP + TDE + covariates
m4_list <- list() # Model 4: Model 3 + district fixed effects


# One pass of this loop fits all four models to imputed dataset i.
for (i in 1:5) {
  # Run Models
  # Model 1: Only AP
  # The simplest possible test: is personalism alone related to pork barrel
  # behaviour? Nothing is held constant, so this is a raw association.
  m1_list[[i]] <- fenegbin(lpgbills ~ pers.hat, se = "hetero", data = hond_list[[i]])

  # Model 2: AP + TDE
  # Add interparty competition, to check that AP is not simply standing in for
  # the other dimension of electoral incentives.
  m2_list[[i]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat, se = "hetero", data = hond_list[[i]])

  # Model 3: Add covariates
  # Hold constant the legislator's own characteristics: membership of the
  # governing party, seniority, a seat on the chamber's board, a committee
  # chair, membership of a small party, and sex.
  m3_list[[i]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[i]])

  # Model 4: Add district fixed effects using fixed-effects negative binomial
  # The strictest specification: everything in Model 3, plus a separate
  # intercept per district (the terms after the | bar).
  m4_list[[i]] <- fenegbin(lpgbills ~ pers.hat + totalEff.hat + govdeps + seniority + boardm + permcomm + smallparties + sex | district, se = "hetero", data = hond_list[[i]])
}

# Pool estimates using Rubin’s rules across the 5 imputed datasets
# Use mice::pool to combine results
# pool() takes the five fitted models and returns one combined set of
# coefficients and standard errors, applying the same two-part variance
# formula spelled out in the bootstrap function above. summary() then presents
# the pooled result as a tidy table with estimates, standard errors, and
# p-values -- the input `extract_Estimates()` expects.
out1 <- summary(mice::pool(m1_list))
out2 <- summary(mice::pool(m2_list))
out3 <- summary(mice::pool(m3_list))
out4 <- summary(mice::pool(m4_list))

## Table 14.1: Association between AP and Initiated Local Goods Bills

# Collect pooled estimates for each model
# Put the four pooled summaries in a list, then run extract_Estimates() on each
# one. sapply(X, FUN) applies FUN to every element of X; simplify = FALSE keeps
# the four results as separate data frames rather than forcing them together.
model_list <- list(out1, out2, out3, out4)
out_estimates <- sapply(model_list, extract_Estimates, simplify = FALSE)

# Merge results from all models into a single table
# Join model 1 and model 2 on the shared 'variable' column. all = TRUE keeps
# rows present in only one of them -- essential here, because the simpler
# models do not contain the control variables, and those cells must end up
# blank rather than causing the row to disappear.
mytable <- merge(out_estimates[[1]], out_estimates[[2]], by = "variable", all = TRUE)
# Then fold in models 3 and 4 one at a time. This one-line for loop has no
# braces, which is legal when the body is a single statement.
for (i in c(3:4)) mytable <- merge(mytable, out_estimates[[i]], by = "variable", all = TRUE)

# Reorder the rows so that variables appear in a logical sequence
# merge() sorts rows alphabetically by variable name, which is not the order
# readers expect. This line lists the row positions explicitly, putting AP and
# TDE at the top and the constant and dispersion parameter at the bottom.
mytable <- mytable[c(6, 10, 4, 7, 3, 5, 9, 8, 2, 1), ]

# Rename rows using variable names as row names
# Move the variable names into the row names, then delete the now-redundant
# column. Assigning NULL to a column is how R removes it.
row.names(mytable) <- mytable$variable
mytable$variable <- NULL

# Calculate Adjusted Pseudo R2 for each model (average across imputations)
# Two goodness-of-fit measures are reported beneath the coefficients.
# Reading this line from the inside out:
#   lapply(m1_list, function(x) x$pseudo_r2)  pulls the pseudo-R2 out of each
#     of the five fitted models, returning a list of five numbers;
#   Reduce('+', ...)  adds that list up into a single total;
#   / length(m1_list) divides by five to give the average;
#   round(..., 3)     keeps three decimal places.
r2_m1 <- round(Reduce("+", lapply(m1_list, function(x) x$pseudo_r2)) / length(m1_list), 3)
r2_m2 <- round(Reduce("+", lapply(m2_list, function(x) x$pseudo_r2)) / length(m2_list), 3)
r2_m3 <- round(Reduce("+", lapply(m3_list, function(x) x$pseudo_r2)) / length(m3_list), 3)
r2_m4 <- round(Reduce("+", lapply(m4_list, function(x) x$pseudo_r2)) / length(m4_list), 3)

# Calculate Bayesian Information Criterion (BIC) for each model
# BIC rewards fit and penalises complexity: -2 x log-likelihood, plus a penalty
# of log(number of observations) for each parameter estimated. Lower is better.
# It is computed by hand here because the quantity has to be averaged across
# the five imputations in the same way as the pseudo-R2 above.
bic_m1 <- round(Reduce("+", lapply(m1_list, function(x) -2 * x$loglik + log(x$nobs) * x$nparams)) / length(m1_list), 3)
bic_m2 <- round(Reduce("+", lapply(m2_list, function(x) -2 * x$loglik + log(x$nobs) * x$nparams)) / length(m2_list), 3)
bic_m3 <- round(Reduce("+", lapply(m3_list, function(x) -2 * x$loglik + log(x$nobs) * x$nparams)) / length(m3_list), 3)
bic_m4 <- round(Reduce("+", lapply(m4_list, function(x) -2 * x$loglik + log(x$nobs) * x$nparams)) / length(m4_list), 3)

# Add model diagnostics and metadata as rows in the table
# rbind() appends a row to the bottom of the table; the text before the = sign
# becomes that row's name. Four values, one per model column.
mytable <- rbind(mytable, "Fixed Effects by District" = c("No", "No", "No", "Yes"))
mytable <- rbind(mytable, "Adj. Pseudo R2" = c(r2_m1, r2_m2, r2_m3, r2_m4))
mytable <- rbind(mytable, "BIC" = c(bic_m1, bic_m2, bic_m3, bic_m4))
mytable <- rbind(mytable, "Observations" = c(m1_list[[1]]$nobs, m2_list[[1]]$nobs, m3_list[[1]]$nobs, m4_list[[1]]$nobs))

# Rename rows to match presentation in book
# Swap the raw variable names (pers.hat, govdeps, ...) for the labels used in
# the printed book. The order here must match the row order set above.
# 'Theta' is the negative binomial dispersion parameter: it measures how much
# more variable the bill counts are than a Poisson model would allow.
row.names(mytable) <- c(
  "AP",
  "TDE",
  "Government's party member",
  "Seniority",
  "Board member",
  "Committee chair",
  "Small Party",
  "Sex = Female",
  "Constant",
  "Theta",
  "Fixed Effects by District",
  "Adj. Pseudo R2",
  "BIC",
  "Observations"
)

# Replace NAs with empty strings
# NA is R's marker for a missing value. Here the NAs are structural -- a
# control variable simply does not appear in Models 1 and 2 -- so they are
# blanked out rather than printed. is.na() returns TRUE wherever a cell is
# missing, and those cells are overwritten with the empty string.
mytable[is.na(mytable)] <- ""

# Create table
# kbl() is kableExtra's table renderer. caption sets the title, col.names
# relabels the four columns as (1) to (4) in the usual journal style, and
# format = 'html' produces the HTML seen on this page.
mytable %>%
  kbl(
    caption = "Association between AP and the Number of Initiated Local Goods Bills",
    col.names = c("(1)", "(2)", "(3)", "(4)"),
    format = "html",
    linesep = "",
  )

## Prepare data: average AP and TDE by district magnitude and ballot type
# Extract unique rows for relevant variables from one imputed dataset
# Create a small dataset and split by system (CLPR and FrLPR)

# The analysis data have one row per legislator-year, but AP and TDE vary only
# by system and district magnitude. unique() collapses the duplicates, leaving
# 22 distinct rows: eleven magnitudes under each of the two ballot structures.
smallData <- unique(hond_list[[1]][, c("pers.hat", "M", "tballot", "totalEff.hat")])
smallData <- smallData[order(smallData$tballot, smallData$pers.hat), ]

# Separate datasets for pre-reform (CLPR) and post-reform (FR-LPR)
# CLPR  = closed-list PR, the pre-2004 system  (tballot == 0)
# FR-LPR = free-list PR, the post-2004 system  (tballot == 1)
clpr <- unique(smallData[smallData$tballot == 0, ])
flpr <- smallData[smallData$tballot == 1, ]

# Add correct AP and TDE (average across all five datasets)
# Delete the single-dataset predictions (assigning NULL removes a column; the
# chained <- removes both in one statement), then merge in the values averaged
# across all five imputations, which were computed earlier as avg_TDE_AP.
smallData$totalEff.hat <- smallData$pers.hat <- NULL
smallData <- merge(smallData, avg_TDE_AP, by.x = c("M", "tballot"), by.y = c("M", "tballot"))

# Split data
# Four plain numeric vectors, each sorted by district magnitude so that the
# entries line up position by position: clpr_ap[3] and flpr_ap[3] describe the
# same district size before and after the reform.
clpr_ap <- clpr$pers.hat[order(clpr$M)]
flpr_ap <- flpr$pers.hat[order(flpr$M)]
clpr_tde <- clpr$totalEff.hat[order(clpr$M)]
flpr_tde <- flpr$totalEff.hat[order(flpr$M)]


### Panel (a): Average Personalism by district magnitude and ballot type ###
# This commented-out line would open a PDF device and write the figure to a
# file instead of showing it on screen. It pairs with the dev.off() at the end
# of the block. Both are disabled so the plot appears inline in this page.
# pdf('~/Downloads/Plots/APHonduras.pdf', width = 10)


# Build dataset for plotting
# Stack the pre- and post-reform AP values into one column, with the matching
# magnitudes beside them. cbind() ("column bind") glues vectors together side
# by side; sort(clpr$M) is repeated because the same eleven magnitudes apply to
# both eras.
ap <- as.data.frame(cbind(AP = c(clpr_ap, flpr_ap), M = c(sort(clpr$M), sort(clpr$M))))
# cbind() can turn numbers into text, so both columns are converted back to
# numeric to be safe.
ap$M <- as.numeric(ap$M)
ap$AP <- as.numeric(ap$AP)
# Sorting by magnitude interleaves the two eras, so the bars come out in
# pre/post pairs: M=1 before, M=1 after, M=2 before, M=2 after, and so on.
ap <- ap[order(ap$M), ]

# Plot AP barplot
# Margins: bottom 6, left 7 (room for the long y-axis label), top 5, right 0.
par(mar = c(6, 7, 5, 0))
# barplot() draws the bars and, usefully, RETURNS the horizontal position of
# each bar. Storing that in difAPplot is what makes it possible to place the
# axis labels and dividing lines at exactly the right coordinates below.
difAPplot <- barplot(ap$AP - 0.3, # Offset to align baseline with y=0.3
  col = c("grey70", "grey50"), border = FALSE, # Grey70 = CLPR, Grey50 = FR-LPR
  names.arg = "", ylim = c(0, 0.25), cex.axis = 1.5,
  ylab = c("Average Personalism"), cex.lab = 2,
  cex.names = 1.5, xlab = "District Magnitude", axes = FALSE
)
# Notes on those arguments:
#   ap$AP - 0.3  subtracts 0.3 from every bar height. All AP values exceed 0.3,
#                so full-height bars would differ only in their tips; cutting
#                the axis at 0.3 makes the differences visible. The true values
#                are restored on the axis labels immediately below.
#   col = c(...) recycles two greys across the bars, which is why the pre/post
#                interleaving arranged above produces a light/dark pairing.
#   names.arg = '' suppresses automatic labels; they are added by hand.
#   cex.axis / cex.lab / cex.names are text magnification factors.
#   axes = FALSE turns off the default axes so custom ones can be drawn.

# Customize y-axis (adjusted for 0.3 offset)
# axis(2, ...) draws the LEFT axis (1 = bottom, 2 = left, 3 = top, 4 = right).
# The tick marks sit at the plotted heights, but the printed labels have the
# 0.3 added back, so the reader sees true AP values.
axis(2,
  at = c(0, 0.05, 0.10, 0.15, 0.2, 0.25),
  labels = c(0, 0.05, 0.10, 0.15, 0.2, 0.25) + 0.3,
  cex.axis = 1.5
)

# Add district magnitude labels to x-axis (every other tick for clarity)
# seq(1, 22, by = 2) gives 1, 3, 5, ... 21: the first bar of each pair. Adding
# 0.6 nudges the label to sit between the two bars of the pair.
# tick = FALSE prints the labels without drawing tick marks.
axis(1,
  at = difAPplot[seq(1, 22, by = 2), 1] + 0.6,
  labels = sort(clpr$M), cex.axis = 1.65, tick = FALSE
)

# Draw vertical dotted lines between grouped bars (CLPR-FRPR)
# Compute where each divider goes, then draw them one at a time.
# segments() draws a straight line from (x0, y0) to (x1, y1); holding x0 equal
# to x1 makes it vertical, and lty = 2 makes it dashed.
forLines <- difAPplot[seq(2, 21, by = 2), 1] + 0.6
for (i in 1:length(forLines)) segments(x0 = forLines[i], x1 = forLines[i], y0 = 0, y1 = 0.3, lty = 2)

# Add bottom axis tick marks (blank labels for spacing)
# A second bottom axis, this time with tick marks but no text -- rep('', 12)
# repeats the empty string twelve times -- so each pair of bars is visually
# bracketed.
axis(1,
  at = c(0, forLines, 25.9 + 0.6),
  labels = rep("", 12), cex.axis = 2, tick = TRUE
)

# Add legend
# new = TRUE tells R to draw over the existing plot rather than start a fresh
# one, which lets the legend be positioned in the top margin.
par(mar = c(0, 0, 0, 0), new = TRUE)
# bty = 'n' suppresses the box around the legend; border = F removes the
# outlines of the colour swatches.
legend(ncol = 2, x = 5, y = .32, legend = c("Pre-Reform", "Post-Reform"), fill = c("grey70", "grey50"), border = F, bty = "n", cex = 1.55)

# Closes the PDF device opened by the pdf() line above; disabled here for the
# same reason.
# dev.off()

### Panel (b): Total Duvergerian Effect ###
# pdf('~/Downloads/Plots/TDEHonduras.pdf', width = 10)

# Same construction as panel (a): stack pre- and post-reform values, then sort
# by magnitude so the bars interleave into light/dark pairs.
tde <- as.data.frame(cbind(TDE = c(clpr_tde, flpr_tde), M = c(sort(clpr$M), sort(clpr$M))))
tde <- tde[order(tde$M), ]
par(mar = c(6, 7, 5, 0))
# No offset is subtracted here, and axes = T keeps the default vertical axis,
# because TDE values start near zero and need no truncation to be legible.
difTDEplot <- barplot(tde$TDE,
  col = c("grey70", "grey50"), border = FALSE,
  names.arg = "", ylim = c(0, 3.5), cex.axis = 1.5,
  ylab = c("Total Duvergerian Effect"), cex.lab = 2,
  cex.names = 1.5, xlab = "District Magnitude", axes = T
)

# Add magnitude labels on bottom axis
# Again positioned from the bar coordinates returned by barplot().
axis(1,
  at = difTDEplot[seq(1, 22, by = 2), 1] + 0.6,
  labels = sort(clpr$M), cex.axis = 1.65, tick = FALSE
)

# Vertical dotted lines
# Drawn up to 3.5 here, matching this panel's vertical range.
forLines <- difTDEplot[seq(2, 21, by = 2), 1] + 0.6
for (i in 1:length(forLines)) segments(x0 = forLines[i], x1 = forLines[i], y0 = 0, y1 = 3.5, lty = 2)

# Add bottom ticks
axis(1,
  at = c(0, forLines, 25.9 + 0.6),
  labels = rep("", 12), cex.axis = 2, tick = TRUE
)
par(mar = c(0, 0, 0, 0), new = TRUE)

# Add legend
legend(ncol = 2, x = 5, y = 4.5, legend = c("Pre-Reform", "Post-Reform"), fill = c("grey70", "grey50"), border = F, bty = "n", cex = 1.55)

# dev.off()

## Calculate Expected Values
expValues <- list()
# One pass per district magnitude. clpr_ap[i] and flpr_ap[i] are the pre- and
# post-reform AP values for the same district size, which is exactly the pair
# boot.predict.pork() needs. length(clpr_ap) is used instead of a hard-coded 11
# so the loop still works if the set of magnitudes ever changes.
for (i in 1:length(clpr_ap)) {
  # For each district magnitude, simulate expected bills under CLPR and FR-LPR
  expValues[[i]] <- boot.predict.pork(myvalueCLPR = clpr_ap[i], myvalueFLPR = flpr_ap[i], reps = 500)
}

# Convert list output into a data frame
# Each element of expValues is a six-number vector. rbind() stacks them into a
# table with one row per magnitude and six columns, in the order the function
# named them: CLPR estimate, CLPR error, FLPR estimate, FLPR error, difference,
# difference error.
expValuesNew <- as.data.frame(do.call("rbind", expValues))

# Compute 95% confidence intervals for the expected difference in counts
# Column 5 is the estimated difference and column 6 its standard error.
# qnorm(0.975) is 1.96, the multiplier that puts 95% of a normal distribution
# inside the interval. ub = upper bound, lb = lower bound; these become the
# error bars drawn on panel (a) below.
expValuesNew$ub <- expValuesNew[, 5] + qnorm(0.975) * expValuesNew[, 6]
expValuesNew$lb <- expValuesNew[, 5] - qnorm(0.975) * expValuesNew[, 6]

# For interpretation used in the chapter: Effect size interpretation in standardized terms

# Dividing an effect by the standard deviation of the outcome expresses it in
# "standard deviation units", the common currency of effect sizes. sd() computes
# the standard deviation of the observed bill counts.
# Note that rows are addressed by number and columns by the names assigned
# inside boot.predict.pork(). Row 10 and row 3 correspond to M = 20 and M = 3
# because expValuesNew was built in ascending order of district magnitude.

# Calculate standardized effect size for M = 20
expValuesNew[10, "Point Estimate - Diff"] / sd(hond$lpgbills)

# Calculate standardized effect size for M = 3
expValuesNew[3, "Point Estimate - Diff"] / sd(hond$lpgbills)

# pdf('~/Downloads/PorkBarrel_CareyShugart.pdf')
# No Controls

# Fit a negative binomial model with an interaction between M and ballot_type
# Note that this uses the observed ballot type and district magnitude directly,
# rather than the GBM-predicted AP used everywhere else -- the point is to test
# Carey and Shugart's argument on its own terms.
c_and_s_noControls <- fenegbin(lpgbills ~ M * ballot_type, se = "hetero", data = hond_list[[1]])
# Print summary of the model, including coefficients, standard errors, and significance
# The row to watch is the interaction, printed as M:ballot_type.
summary(c_and_s_noControls)

# Plot the interaction effect of M and ballot_type using sjPlot
# type = "int" asks plot_model() for an interaction plot: predicted outcome on
# the vertical axis, the first term in `terms` along the horizontal axis, and a
# separate line for each value of the second term. Parallel lines mean no
# interaction; converging or diverging lines mean the effect of one variable
# depends on the other.
# robust = TRUE uses the heteroskedasticity-robust standard errors when drawing
# the confidence bands.
sjPlot::plot_model(c_and_s_noControls,
  type = "int", terms = c("M", "ballot_type"), robust = TRUE,
  axis.title = c("M", "Attention to Candidate"),
  title = "Expected Value of Attention to Candidate - No Controls",
  legend.title = "Ballot Type"
)

# With Controls

# Fit a negative binomial model with interaction + control variables
# Identical to the previous model except for the six legislator-level controls
# added after the interaction. Keeping everything else the same is what makes
# the two sets of results directly comparable.
c_and_s_wControls <- fenegbin(lpgbills ~ M * ballot_type + govdeps + seniority + boardm + permcomm + smallparties + sex, se = "hetero", data = hond_list[[1]])

# Show coefficient estimates, robust SEs, and significance
summary(c_and_s_wControls)

# Plot interaction effect again, this time accounting for controls
# When controls are present, plot_model() holds them at typical values while
# varying M and ballot type -- the same "hold everything else constant" logic
# used inside boot.predict.pork().
sjPlot::plot_model(c_and_s_wControls,
  type = "int", terms = c("M", "ballot_type"), robust = TRUE,
  axis.title = c("M", "Attention to Candidate"),
  title = "Expected Value of Attention to Candidate - With Controls",
  legend.title = "Ballot Type"
)

# dev.off()