# Load relevant libraries

library(tidyverse) # bundle of data-handling packages (dplyr, ggplot2, ...)
library(haven) # reads Stata .dta files, preserving value labels
library(MASS) # statistical functions; NOTE it masks dplyr::select()
library(dplyr) # re-attached here, but see the note above on masking
library(sjPlot) # model plots and tables
library(lme4) # multilevel / mixed-effects models
library(gbm) # gradient boosting machines: supply the AP and TDE predictions
library(gridExtra) # arranges several plot panels into one figure
library(stargazer) # formats regression results into publication tables

# Set the path to Chapter 12 materials

# the short file names below. See the note above on why this is repeated.


# The Van Dusky-Allen & Touchton dataset, in Stata format. Supplies the
# chamber-level variables, including numscl -- the number of standing
# committees, which is the outcome for the first half of the chapter.
dta <- read_dta("data/ch12/dta/VanDuskyAllenTouchtonPRQ.dta")

# National GBM
# ---------------------------------------------------------------------------
# Pre-computed country-level TDE and AP scores.
#
# THIS CHAPTER DELIBERATELY USES A DIFFERENT SCORES FILE FROM CHAPTERS 6 TO 11.
#
# Those chapters read RealSystems_Scores_GBM_Aug_2024.RData. This one reads the
# earlier RealSystems_Scores_GBM.RData, and must continue to, because it is the
# file the published Table 12.1 was computed from.
#
# The distinction matters here and not elsewhere because this chapter uses the
# stored scores DIRECTLY as its explanatory variables -- see the merge below,
# where pers.hat1 becomes AP and totalEff.hat1 becomes TDE. The two files carry
# genuinely different scores, so switching moves every coefficient in
# Table 12.1. Tested: the August 2024 file gives AP 1.550 / TDE 0.321 in
# model 1 against the published 1.556 / 0.290. Substantive conclusions are
# unaffected either way, but the printed table would no longer match.
#
# Chapter 11 could be switched safely precisely because it does NOT use the
# stored scores -- it recomputes them from the GBMs.
# ---------------------------------------------------------------------------
load("data/shared/RealSystems_Scores_GBM.RData")

# District GBM
# The fitted models from Chapter 6, used in the second half of the chapter to
# score individual districts rather than whole countries.
load(file = "data/shared/AP_district_objects_t515_d7.RData")
load(file = "data/shared/TDE_district_objects_t305_d15.RData")

# Work with a smaller set of variables from VDAT
dta_small <- dta %>%
  dplyr::select(
    year, country_name,
    numscl, v2psprlnks_osp,
    numberlower, eleclower,
    executive, effn, bicameral,
    restrictionappt, logg
  ) %>%
  rename(country = country_name)

dta_small$country <- dplyr::recode(dta_small$country,
  "Korea South" = "South Korea",
  "Slovak Republic" = "Slovakia",
  "United States" = "United States of America"
)

# Merging the two datasets

# ---------------------------------------------------------------------------
# Keep only the columns this chapter needs from the scores file, BEFORE merging.
#
# This is a safeguard rather than a fix for a present problem. drop_na() a few
# lines below removes any row with a missing value in ANY column of the merged
# table -- including columns that play no part in the analysis. The file loaded
# above has no gaps in the columns kept here, so this line changes nothing
# today: the sample is 70 chambers with or without it.
#
# It matters because that guarantee is fragile. A later scores file with an
# extra column carrying missing values would silently delete observations. The
# August 2024 file used by chapters 6 to 11 is exactly such a file -- it adds
# ten columns, three with gaps (`democ`, derived from Polity, has 37), which
# would drop rows on the strength of variables this chapter never reads.
#
# The 14 columns kept are the merge key, the year used for the date filter, the
# five imputed TDE and AP predictions, and their averages. Nothing else in this
# chapter touches the scores file.
# ---------------------------------------------------------------------------
scoreCols <- c(
  "country", "year",
  paste0("totalEff.hat", 1:5),
  paste0("pers.hat", 1:5),
  "totalEff.hat", "pers.hat"
)

dta_merged <- merge(data2export[, scoreCols], dta_small,
  by = "country",
  all.y = T
)

# Remove pre-2005, post-2015 elections
# `year.x` is the election year from the scores file; `year.y` is VDAT's. Only
# `country` and `year` appear in both inputs, so those are the only names the
# merge disambiguates with .x / .y suffixes.
dta_merged <- dta_merged[dta_merged$year.x >= 2005 & dta_merged$year.x <= 2015, ]

# Drop rows with any missing value. See the note above on why the column
# selection has to happen before this line rather than after it.
dta_merged <- dta_merged %>% drop_na()

dta_merged <- dta_merged %>% arrange(country, year.x)

# We want to keep from our dataset the election immediately prior to 2015
# (2015 is the year for which VDAT have data)
# The following snippet drops information from all elections
# except the one immediately prior to 2015
dta_merged$proper_elections <- 0
for (i in 1:(nrow(dta_merged) - 1)) {
  if (dta_merged$country[i] == dta_merged$country[i + 1]) {
    dta_merged$proper_elections[i] <- 1
  }
}
dta_merged <- dta_merged %>%
  dplyr::filter(proper_elections != 1) %>%
  dplyr::select(-proper_elections)

#################
#### Germany ####
#################


# A large majority of MPs (7079) have legislature value 999
# Though about 116 of these have a value as having been elected
# in either a nominal or list district, and even a few appear
# to have some committee assignment, none of these appear to
# be "real" MPs. From now on, we proceed to eliminate them.

ger <- read.csv("data/ch12/csv/germany.csv")

ger.old <- ger
ger <- ger %>%
  dplyr::filter(legislature != 999)

M <- c(unique(ger$magnitude_list), 1)

# Column M guarantees that we are only selecting elected MP
# un-elected candidates have elected_list == 0 and elected_nom == 0
ger <- ger %>%
  mutate(M = if_else(elected_list == 1, magnitude_list, magnitude_nom))

### Get the objects from list

optimalGBMIntra <- totalAP.objects$optimalGBM
optimalGBMInter <- totalEffENP.objects$optimalGBM

# Create ger_real_cases Dataset

ger_real_cases <- rbind(
  data.frame(tballot = 1, ballot_type = "closed", M = M, formula = "saintelague", new.nvotes = "One", pool_level = "party", threshold = 5),
  data.frame(tballot = 0, ballot_type = "closed", M = M, formula = "hare", new.nvotes = "One", pool_level = "party", threshold = 5),
  data.frame(tballot = 2, ballot_type = "closed", M = M, formula = "dhondt", new.nvotes = "One", pool_level = "party", threshold = 5)
)

### Get the objects from list
ger_real_cases$new.nvotes <- factor(ger_real_cases$new.nvotes, levels = optimalGBMIntra[[1]]$var.levels[[3]])
ger_real_cases$pool_level <- factor(ger_real_cases$pool_level, levels = optimalGBMIntra[[1]]$var.levels[[2]])
ger_real_cases$ballot_type <- factor(ger_real_cases$ballot_type, levels = optimalGBMIntra[[1]]$var.levels[[1]])

# Create ger_real_cases Dataset
ger_real_cases$formula <- factor(ger_real_cases$formula, levels = optimalGBMIntra[[1]]$var.levels[[6]])
ger_real_cases$new.nvotes[ger_real_cases$M == 1] <- "One"
ger_real_cases$ballot_type[ger_real_cases$M == 1] <- "closed"
ger_real_cases$formula[ger_real_cases$M == 1] <- "plurality"

# Reorder the data
ger_real_cases <- ger_real_cases[order(ger_real_cases$ballot_type, ger_real_cases$M), ]

# This loop generates five sets of predictions for each MP using both the interparty (totalEff.hat) and intraparty (pers.hat) GBM models. The predictions are stored in the list gerData for further analysis.
gerData <- list()
for (i in 1:5) {
  # Interparty
  ger_real_cases$totalEff.hat <- predict(optimalGBMInter[[i]], ger_real_cases, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  ger_real_cases$pers.hat <- predict(optimalGBMIntra[[i]], ger_real_cases, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  gerData[[i]] <- ger_real_cases
}

ger1 <- gerData[[1]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat1 = totalEff.hat, pers.hat1 = pers.hat)

ger2 <- gerData[[2]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat2 = totalEff.hat, pers.hat2 = pers.hat)

ger3 <- gerData[[3]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat3 = totalEff.hat, pers.hat3 = pers.hat)

ger4 <- gerData[[4]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat4 = totalEff.hat, pers.hat4 = pers.hat)

ger5 <- gerData[[5]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat5 = totalEff.hat, pers.hat5 = pers.hat)

# Combine all the predictions
ger_gbms <- ger1 %>%
  full_join(ger2) %>%
  full_join(ger3) %>%
  full_join(ger4) %>%
  full_join(ger5)

# Merge the predictions with the original data
ger <- merge(ger, ger_gbms, by = "M", all = T)

# The code above generated AP and TDE for every MP under all different
# rules in Germany. The 'drop' column below is to drop all observations
# where a given rule was not applicable in an election-year

ger$drop <- ifelse(ger$formula == "plurality" & ger$tballot == 1, 1,
  ifelse(ger$formula == "dhondt" & ger$election_year == 1983, 1,
    ifelse(ger$formula == "hare" & ger$election_year == 1987, 1,
      ifelse(ger$formula == "hare" & ger$election_year == 1998, 1,
        ifelse(ger$formula == "hare" & ger$election_year == 2005, 1,
          ifelse(ger$formula == "saintelague" & ger$election_year == 2009, 1, 0)
        )
      )
    )
  )
)

ger <- ger %>%
  dplyr::filter(drop == 1, ) %>%
  dplyr::select(-drop)

# Generate the pork variable in accordance to Shugart et al.'s book

ger$pork <- ifelse(ger$sedepe10 > 0 | # sedepe11 is missing
  ger$sedepe18 > 0 |
  ger$sedepe19 > 0 |
  ger$sedepe28 > 0 |
  ger$sedepe36 > 0, 1, 0)
unique(with(ger, rowSums(cbind(pork, pork_committee)))) # This should turn 2 and NA

ger$pork <- replace_na(ger$pork, 0)

# Defining CDU and CSU as a single party

ger$party_name <- ifelse(ger$party_name == "CDU" | ger$party_name == "CSU",
  "CDU/CSU", ger$party_name
)

###############################
#### Japan -- mixed system ####
###############################


jpn <- read.csv("data/ch12/csv/japan.csv")

M <- c(unique(jpn$magnitude_list), 1)

# It appears that JPN has a large number of individuals with elected_list=0 and elected_nom=0
jpn <- jpn %>%
  mutate(M = case_when(
    elected_list == 1 ~ magnitude_list,
    elected_nom == 1 ~ magnitude_nom,
    TRUE ~ NA_integer_
  ))

jpn <- jpn[jpn$M >= 1 & !is.na(jpn$M), ]

M <- M[-1]

# Create jpn_real_cases Dataset

jpn_real_cases <- rbind(data.frame(tballot = 2, ballot_type = "closed", M = M, formula = "dhondt", new.nvotes = "One", pool_level = "party", threshold = 0))
jpn_real_cases$new.nvotes <- factor(jpn_real_cases$new.nvotes, levels = optimalGBMIntra[[1]]$var.levels[[3]])
jpn_real_cases$pool_level <- factor(jpn_real_cases$pool_level, levels = optimalGBMIntra[[1]]$var.levels[[2]])
jpn_real_cases$ballot_type <- factor(jpn_real_cases$ballot_type, levels = optimalGBMIntra[[1]]$var.levels[[1]])
jpn_real_cases$formula <- factor(jpn_real_cases$formula, levels = optimalGBMIntra[[1]]$var.levels[[6]])
jpn_real_cases$new.nvotes[jpn_real_cases$M == 1] <- "One"
jpn_real_cases$ballot_type[jpn_real_cases$M == 1] <- "closed"
jpn_real_cases$formula[jpn_real_cases$M == 1] <- "plurality"

# Reorder the data
jpn_real_cases <- jpn_real_cases[order(jpn_real_cases$ballot_type, jpn_real_cases$M), ]

# In this loop, we generate predictions for interparty and intraparty effects using the pre-trained optimalGBMInter and optimalGBMIntra models. These predictions are stored in jpnData for further analysis.
# Predict Interparty and Intraparty
jpnData <- list()
for (i in 1:5) {
  # Interparty
  jpn_real_cases$totalEff.hat <- predict(optimalGBMInter[[i]], jpn_real_cases, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  jpn_real_cases$pers.hat <- predict(optimalGBMIntra[[i]], jpn_real_cases, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  jpnData[[i]] <- jpn_real_cases
}


# Combine all the predictions
jpn1 <- jpnData[[1]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat1 = totalEff.hat, pers.hat1 = pers.hat)

jpn2 <- jpnData[[2]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat2 = totalEff.hat, pers.hat2 = pers.hat)

jpn3 <- jpnData[[3]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat3 = totalEff.hat, pers.hat3 = pers.hat)

jpn4 <- jpnData[[4]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat4 = totalEff.hat, pers.hat4 = pers.hat)

jpn5 <- jpnData[[5]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat5 = totalEff.hat, pers.hat5 = pers.hat)

# The five datasets (jpn1 through jpn5) are combined into one large dataset
jpn_gbms <- jpn1 %>%
  full_join(jpn2) %>%
  full_join(jpn3) %>%
  full_join(jpn4) %>%
  full_join(jpn5)

# The jpn dataset is merged with the jpn_gbms dataset
jpn <- merge(jpn, jpn_gbms, by = "M", all = T)

# A pork variable is created based on the presence of certain activities (indicated by columns like sedepe10, sedepe11, etc.). If any of these columns have non-zero values, the MP is considered involved in pork-barrel politics.
jpn$pork <- ifelse(jpn$sedepe10 > 0 | # sedepe28 is missing
  jpn$sedepe11 > 0 |
  jpn$sedepe18 > 0 |
  jpn$sedepe19 > 0 |
  jpn$sedepe36 > 0, 1, 0)

jpn$pork <- replace_na(jpn$pork, 0)

# unique (with (jpn, rowSums (cbind (pork, porkcommittee)))) # not possible because porkcommittee is all NA

#######################
#### Japan -- SNTV ####
#######################


jpn1 <- read.csv("data/ch12/csv/japan.csv")

jpn1 <- jpn1[jpn1$election_year < 1996, ] # Only elections before 1996, under the old system

# To make japan sntv a "different" country from japan MM

jpn1$country <- "Japan1"

M <- 1:6

jpn1$M <- jpn1$elected_nom * jpn1$magnitude_nom

jpn1 <- jpn1[jpn1$M >= 1 & !is.na(jpn1$M), ]

# Create jpn1_real_cases Dataset

jpn1_real_cases <- rbind(data.frame(tballot = 2, ballot_type = "open", M = M, formula = "plurality", new.nvotes = "One", pool_level = "candidate", threshold = 0))
jpn1_real_cases$new.nvotes <- factor(jpn1_real_cases$new.nvotes, levels = optimalGBMIntra[[1]]$var.levels[[3]])
jpn1_real_cases$pool_level <- factor(jpn1_real_cases$pool_level, levels = optimalGBMIntra[[1]]$var.levels[[2]])
jpn1_real_cases$ballot_type <- factor(jpn1_real_cases$ballot_type, levels = optimalGBMIntra[[1]]$var.levels[[1]])
jpn1_real_cases$formula <- factor(jpn1_real_cases$formula, levels = optimalGBMIntra[[1]]$var.levels[[6]])

# Reorder the data
jpn1_real_cases <- jpn1_real_cases[order(jpn1_real_cases$ballot_type, jpn1_real_cases$M), ]

### We need to run predictions five times
# Predict Interparty and Intraparty
jpn1Data <- list()
for (i in 1:5) {
  # Interparty
  jpn1_real_cases$totalEff.hat <- predict(optimalGBMInter[[i]], jpn1_real_cases, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  jpn1_real_cases$pers.hat <- predict(optimalGBMIntra[[i]], jpn1_real_cases, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  jpn1Data[[i]] <- jpn1_real_cases
}

jpn11 <- jpn1Data[[1]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat1 = totalEff.hat, pers.hat1 = pers.hat)

jpn12 <- jpn1Data[[2]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat2 = totalEff.hat, pers.hat2 = pers.hat)

jpn13 <- jpn1Data[[3]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat3 = totalEff.hat, pers.hat3 = pers.hat)

jpn14 <- jpn1Data[[4]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat4 = totalEff.hat, pers.hat4 = pers.hat)

jpn15 <- jpn1Data[[5]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat5 = totalEff.hat, pers.hat5 = pers.hat)

# Combine all the predictions
jpn1_gbms <- jpn11 %>%
  full_join(jpn12) %>%
  full_join(jpn13) %>%
  full_join(jpn14) %>%
  full_join(jpn15)

jpn1 <- merge(jpn1, jpn1_gbms, by = "M", all = T)

# The pork variable is created in the same way as for the mixed system, identifying MPs involved in pork-barrel activities.
jpn1$pork <- ifelse(jpn1$sedepe10 > 0 | # sedepe28 missing
  jpn1$sedepe11 > 0 |
  jpn1$sedepe18 > 0 |
  jpn1$sedepe19 > 0 |
  jpn1$sedepe36 > 0, 1, 0)

jpn1$pork <- replace_na(jpn1$pork, 0)

# unique (with (jpn1, rowSums (cbind (pork, porkcommittee))))

#####################
#### New Zealand ####
#####################


nzl <- read.csv("data/ch12/csv/new_zealand.csv")

# M is created to represent the electoral magnitude based on whether a candidate was elected through a list or nominal system. Non-elected candidates are excluded from the dataset.
M <- c(1, 50, 51, 53, 55)

nzl <- nzl %>%
  mutate(M = case_when(
    elected_list == 1 ~ magnitude_list,
    elected_nom == 1 ~ magnitude_nom,
    TRUE ~ NA_integer_
  ))

nzl <- nzl[nzl$M >= 1 & !is.na(nzl$M), ]

# nzl <- nzl[nzl$election_year==1996 | nzl$election_year==1999 | nzl$election_year==2002 | nzl$election_year==2005 | nzl$election_year==2008,]


# This line ensures that only election years from 1996 onward are included, reflecting the adoption of the MMP system in New Zealand starting in 1996.

nzl <- nzl[nzl$election_year >= 1996, ]

# Create nzl_real_cases Dataset

nzl_real_cases <- rbind(data.frame(data.frame(tballot = 0, ballot_type = "closed", M = M, formula = "saintelague", new.nvotes = "One", pool_level = "party", threshold = 5)))
nzl_real_cases$new.nvotes <- factor(nzl_real_cases$new.nvotes, levels = optimalGBMIntra[[1]]$var.levels[[3]])
nzl_real_cases$pool_level <- factor(nzl_real_cases$pool_level, levels = optimalGBMIntra[[1]]$var.levels[[2]])
nzl_real_cases$ballot_type <- factor(nzl_real_cases$ballot_type, levels = optimalGBMIntra[[1]]$var.levels[[1]])
nzl_real_cases$formula <- factor(nzl_real_cases$formula, levels = optimalGBMIntra[[1]]$var.levels[[6]])
nzl_real_cases$new.nvotes[nzl_real_cases$M == 1] <- "One"
nzl_real_cases$ballot_type[nzl_real_cases$M == 1] <- "closed"
nzl_real_cases$formula[nzl_real_cases$M == 1] <- "plurality"

# The data is ordered by ballot_type and M to facilitate the subsequent modeling process.
nzl_real_cases <- nzl_real_cases[order(nzl_real_cases$ballot_type, nzl_real_cases$M), ]


# In this loop, we generate predictions for both interparty (totalEff.hat) and intraparty (pers.hat) effects using the pre-trained Gradient Boosting Machine (GBM) models.

nzlData <- list()
for (i in 1:5) {
  # Interparty
  nzl_real_cases$totalEff.hat <- predict(optimalGBMInter[[i]], nzl_real_cases, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  nzl_real_cases$pers.hat <- predict(optimalGBMIntra[[i]], nzl_real_cases, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  nzlData[[i]] <- nzl_real_cases
}

# The results from the five iterations in nzlData are processed to select relevant columns and rename the prediction columns for easier identification.

nzl1 <- nzlData[[1]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat1 = totalEff.hat, pers.hat1 = pers.hat)

nzl2 <- nzlData[[2]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat2 = totalEff.hat, pers.hat2 = pers.hat)

nzl3 <- nzlData[[3]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat3 = totalEff.hat, pers.hat3 = pers.hat)

nzl4 <- nzlData[[4]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat4 = totalEff.hat, pers.hat4 = pers.hat)

nzl5 <- nzlData[[5]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat5 = totalEff.hat, pers.hat5 = pers.hat)


# The five datasets (nzl1 through nzl5) are combined into a single dataset, nzl_gbms, which contains the predictions for interparty and intraparty effects for each combination of variables.


nzl_gbms <- nzl1 %>%
  full_join(nzl2) %>%
  full_join(nzl3) %>%
  full_join(nzl4) %>%
  full_join(nzl5)

nzl <- merge(nzl, nzl_gbms, by = "M", all = T)

# The pork variable is created to identify MPs involved in pork-barrel activities.

nzl$pork <- ifelse(nzl$sedepe10 > 0 | # sedepe19, 28, 36
  nzl$sedepe11 > 0 |
  nzl$sedepe18 > 0, 1, 0)

nzl$pork <- replace_na(nzl$pork, 0)


# unique (with (nzl, rowSums (cbind (pork, porkcommittee)))) # here, we have values of 1, which means that we are missing something


#####################
#### New Zealand ####
#####################

# The nzl1 dataset is read, and data for elections prior to 1996 is selected, as this reflects the SNTV system before the switch to MMP.

nzl1 <- read.csv("data/ch12/csv/new_zealand.csv")

nzl1 <- nzl1[nzl1$election_year < 1996, ]

nzl1$country <- "New Zealand1"


# A new variable M is created, reflecting the electoral magnitude based on whether candidates were elected through a list or nominal system.

M <- 1

nzl1 <- nzl1 %>%
  mutate(M = case_when(
    elected_list == 1 ~ magnitude_list,
    elected_nom == 1 ~ magnitude_nom,
    TRUE ~ NA_integer_
  ))

nzl1 <- nzl1[nzl1$M >= 1 & !is.na(nzl1$M), ]


# Create nzl1_real_cases Dataset

nzl1_real_cases <- rbind(data.frame(data.frame(tballot = 0, ballot_type = "closed", M = 1, formula = "plurality", new.nvotes = "One", pool_level = "party", threshold = 0)))
nzl1_real_cases$new.nvotes <- factor(nzl1_real_cases$new.nvotes, levels = optimalGBMIntra[[1]]$var.levels[[3]])
nzl1_real_cases$pool_level <- factor(nzl1_real_cases$pool_level, levels = optimalGBMIntra[[1]]$var.levels[[2]])
nzl1_real_cases$ballot_type <- factor(nzl1_real_cases$ballot_type, levels = optimalGBMIntra[[1]]$var.levels[[1]])
nzl1_real_cases$formula <- factor(nzl1_real_cases$formula, levels = optimalGBMIntra[[1]]$var.levels[[6]])
nzl1_real_cases$new.nvotes[nzl1_real_cases$M == 1] <- "One"
nzl1_real_cases$ballot_type[nzl1_real_cases$M == 1] <- "closed"
nzl1_real_cases$formula[nzl1_real_cases$M == 1] <- "plurality"

# Reorder the data
nzl1_real_cases <- nzl1_real_cases[order(nzl1_real_cases$ballot_type, nzl1_real_cases$M), ]

# Predictions for interparty and intraparty effects are made using the pre-trained models for the SNTV system

nzl1Data <- list()
for (i in 1:5) {
  # Interparty
  nzl1_real_cases$totalEff.hat <- predict(optimalGBMInter[[i]], nzl1_real_cases, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  nzl1_real_cases$pers.hat <- predict(optimalGBMIntra[[i]], nzl1_real_cases, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  nzl1Data[[i]] <- nzl1_real_cases
}

nzl11 <- nzl1Data[[1]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat1 = totalEff.hat, pers.hat1 = pers.hat)

nzl12 <- nzl1Data[[2]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat2 = totalEff.hat, pers.hat2 = pers.hat)

nzl13 <- nzl1Data[[3]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat3 = totalEff.hat, pers.hat3 = pers.hat)

nzl14 <- nzl1Data[[4]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat4 = totalEff.hat, pers.hat4 = pers.hat)

nzl15 <- nzl1Data[[5]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat5 = totalEff.hat, pers.hat5 = pers.hat)


# Combine all the predictions
nzl1_gbms <- nzl11 %>%
  full_join(nzl12) %>%
  full_join(nzl13) %>%
  full_join(nzl14) %>%
  full_join(nzl15)

nzl1 <- merge(nzl1, nzl1_gbms, by = "M", all = T)

# The results from all five iterations are combined and merged with the original dataset
# The pork variable is created for the SNTV system to track MPs involved in pork-barrel activities.
nzl1$pork <- ifelse(nzl1$sedepe10 > 0 |
  nzl1$sedepe11 > 0 |
  nzl1$sedepe18 > 0, 1, 0)

nzl1$pork <- replace_na(nzl1$pork, 0)

# unique (with (nzl1, rowSums (cbind (pork, porkcommittee)))) # here, we have values of 1, which means that we are missing something

#################
#### Bolivia ####
#################

bol <- read.csv("data/ch12/csv/bolivia.csv")

# Create the M variable based on magnitude list
M <- c(unique(bol$magnitude_list), 1)

bol <- bol %>%
  mutate(M = case_when(
    elected_list == 1 ~ magnitude_list,
    elected_nom == 1 ~ magnitude_nom,
    TRUE ~ NA_integer_
  ))

bol <- bol[bol$M >= 1 & !is.na(bol$M), ] # Remove non-elected candidates

# Filter for election years from 1997 onwards
bol <- bol[bol$electionyear >= 1997, ]

# Create the bol_real_cases dataset for various ballot types and formulas

# Set factor levels based on pre-defined model levels
bol_real_cases <- rbind(data.frame(data.frame(tballot = 0, ballot_type = "closed", M = M, formula = "dhondt", new.nvotes = "One", pool_level = "party", threshold = 3)))
bol_real_cases$new.nvotes <- factor(bol_real_cases$new.nvotes, levels = optimalGBMIntra[[1]]$var.levels[[3]])
bol_real_cases$pool_level <- factor(bol_real_cases$pool_level, levels = optimalGBMIntra[[1]]$var.levels[[2]])
bol_real_cases$ballot_type <- factor(bol_real_cases$ballot_type, levels = optimalGBMIntra[[1]]$var.levels[[1]])

# Adjust values for M==1
bol_real_cases$formula <- factor(bol_real_cases$formula, levels = optimalGBMIntra[[1]]$var.levels[[6]])
bol_real_cases$new.nvotes[bol_real_cases$M == 1] <- "One"
bol_real_cases$ballot_type[bol_real_cases$M == 1] <- "closed"
bol_real_cases$formula[bol_real_cases$M == 1] <- "plurality"

# Reorder the data
bol_real_cases <- bol_real_cases[order(bol_real_cases$ballot_type, bol_real_cases$M), ]

# Loop to predict interparty and intraparty effects five times
bolData <- list()
for (i in 1:5) {
  # Interparty
  bol_real_cases$totalEff.hat <- predict(optimalGBMInter[[i]], bol_real_cases, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  bol_real_cases$pers.hat <- predict(optimalGBMIntra[[i]], bol_real_cases, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  bolData[[i]] <- bol_real_cases
}


# Select relevant columns for each prediction iteration and rename
bol1 <- bolData[[1]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat1 = totalEff.hat, pers.hat1 = pers.hat)

bol2 <- bolData[[2]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat2 = totalEff.hat, pers.hat2 = pers.hat)

bol3 <- bolData[[3]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat3 = totalEff.hat, pers.hat3 = pers.hat)

bol4 <- bolData[[4]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat4 = totalEff.hat, pers.hat4 = pers.hat)

bol5 <- bolData[[5]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat5 = totalEff.hat, pers.hat5 = pers.hat)

# Combine all predictions into a single dataset
bol_gbms <- bol1 %>%
  full_join(bol2) %>%
  full_join(bol3) %>%
  full_join(bol4) %>%
  full_join(bol5)

# Merge the predictions with the original dataset
bol <- merge(bol, bol_gbms, by = "M", all = T)

# Create the pork variable based on specific columns indicating pork-barrel politics
bol$pork <- ifelse(bol$sedepe10 > 0 | # sedepe11, 19, 36 missing
  bol$sedepe28 > 0 | # sedepe28 exists, but all 0
  bol$sedepe18 > 0, 1, 0)

bol$pork <- replace_na(bol$pork, 0)

# with (bol, rowSums (cbind (pork, porkcommittee)))

###################
#### Bolivia 1 ####
###################

bol1 <- read.csv("data/ch12/csv/bolivia.csv")

# Filter for specific election years
bol1 <- bol1[bol1$electionyear == 1980 | bol1$electionyear == 1985 | bol1$electionyear == 1989, ]

M <- c(unique(bol1$magnitude_list), 1) # Create M variable for Bolivia1 dataset

# Mark Bolivia1 dataset as a separate country for identification
bol1$country <- "Bolivia1"

# Create M variable for Bolivia1 dataset, based on election method
bol1 <- bol1 %>%
  mutate(M = case_when(
    elected_list == 1 ~ magnitude_list,
    elected_nom == 1 ~ magnitude_nom,
    TRUE ~ NA_integer_
  ))

bol1 <- bol1[bol1$M >= 1 & !is.na(bol1$M), ]


# Create bol_real_cases Dataset

bol1_real_cases <- rbind(data.frame(data.frame(tballot = 0, ballot_type = "closed", M = M, formula = "hare", new.nvotes = "One", pool_level = "party", threshold = 0)))
# Set factor levels based on pre-defined model levels
bol1_real_cases$new.nvotes <- factor(bol1_real_cases$new.nvotes, levels = optimalGBMIntra[[1]]$var.levels[[3]])
bol1_real_cases$pool_level <- factor(bol1_real_cases$pool_level, levels = optimalGBMIntra[[1]]$var.levels[[2]])
bol1_real_cases$ballot_type <- factor(bol1_real_cases$ballot_type, levels = optimalGBMIntra[[1]]$var.levels[[1]])
# Adjust values for M==1
bol1_real_cases$formula <- factor(bol1_real_cases$formula, levels = optimalGBMIntra[[1]]$var.levels[[6]])
bol1_real_cases$new.nvotes[bol1_real_cases$M == 1] <- "One"
bol1_real_cases$ballot_type[bol1_real_cases$M == 1] <- "closed"
bol1_real_cases$formula[bol1_real_cases$M == 1] <- "plurality"

# Reorder the dataset
bol1_real_cases <- bol1_real_cases[order(bol1_real_cases$ballot_type, bol1_real_cases$M), ]

# Loop to predict interparty and intraparty effects five times
bol1Data <- list()
for (i in 1:5) {
  # Interparty
  bol1_real_cases$totalEff.hat <- predict(optimalGBMInter[[i]], bol1_real_cases, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  bol1_real_cases$pers.hat <- predict(optimalGBMIntra[[i]], bol1_real_cases, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  bol1Data[[i]] <- bol1_real_cases
}

# Process the predictions for each iteration and merge the results
bolA <- bol1Data[[1]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat1 = totalEff.hat, pers.hat1 = pers.hat)

bolB <- bol1Data[[2]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat2 = totalEff.hat, pers.hat2 = pers.hat)

bolC <- bol1Data[[3]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat3 = totalEff.hat, pers.hat3 = pers.hat)

bolD <- bol1Data[[4]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat4 = totalEff.hat, pers.hat4 = pers.hat)

bolE <- bol1Data[[5]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat5 = totalEff.hat, pers.hat5 = pers.hat)


# Combine all predictions
bol1_gbms <- bolA %>%
  full_join(bolB) %>%
  full_join(bolC) %>%
  full_join(bolD) %>%
  full_join(bolE)

# Merge the predictions with the original dataset
bol1 <- merge(bol1, bol1_gbms, by = "M", all = T)

bol1$pork <- ifelse(bol1$sedepe10 > 0 | # sedepe11, 19, 36 missing
  bol1$sedepe28 > 0 |
  bol1$sedepe18 > 0, 1, 0)

# Create the pork variable for Bolivia1
bol1$pork <- replace_na(bol1$pork, 0)

# with (bol1, rowSums (cbind (pork, pork_committee)))

###################
#### Bolivia 2 ####
###################

bol2 <- read.csv("data/ch12/csv/bolivia.csv")

# Filter for the 1993 elections
bol2 <- bol2[bol2$electionyear == 1993, ]

# Create M variable for Bolivia2
M <- c(unique(bol2$magnitude_list), 1)

# Mark Bolivia2 dataset as a separate country for identification
bol2$country <- "Bolivia2"

# Create M variable for Bolivia2 dataset, based on election method
bol2 <- bol2 %>%
  mutate(M = case_when(
    elected_list == 1 ~ magnitude_list,
    elected_nom == 1 ~ magnitude_nom,
    TRUE ~ NA_integer_
  ))

# Filter out non-elected candidates
bol2 <- bol2[bol2$M >= 1 & !is.na(bol2$M), ]


# Create bol_real_cases Dataset

bol2_real_cases <- rbind(data.frame(data.frame(tballot = 0, ballot_type = "closed", M = M, formula = "saintelague", new.nvotes = "One", pool_level = "party", threshold = 0)))
# Set factor levels based on pre-defined model levels
bol2_real_cases$new.nvotes <- factor(bol2_real_cases$new.nvotes, levels = optimalGBMIntra[[1]]$var.levels[[3]])
bol2_real_cases$pool_level <- factor(bol2_real_cases$pool_level, levels = optimalGBMIntra[[1]]$var.levels[[2]])
bol2_real_cases$ballot_type <- factor(bol2_real_cases$ballot_type, levels = optimalGBMIntra[[1]]$var.levels[[1]])

# Adjust values for M==1
bol2_real_cases$formula <- factor(bol2_real_cases$formula, levels = optimalGBMIntra[[1]]$var.levels[[6]])
bol2_real_cases$new.nvotes[bol2_real_cases$M == 1] <- "One"
bol2_real_cases$ballot_type[bol2_real_cases$M == 1] <- "closed"
bol2_real_cases$formula[bol2_real_cases$M == 1] <- "plurality"

# Reorder the data
bol2_real_cases <- bol2_real_cases[order(bol2_real_cases$ballot_type, bol2_real_cases$M), ]


# Loop to predict interparty and intraparty effects five times
bol2Data <- list()
for (i in 1:5) {
  # Interparty
  bol2_real_cases$totalEff.hat <- predict(optimalGBMInter[[i]], bol2_real_cases, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  bol2_real_cases$pers.hat <- predict(optimalGBMIntra[[i]], bol2_real_cases, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  bol2Data[[i]] <- bol2_real_cases
}


# Process the predictions for each iteration and merge the results
bola <- bol2Data[[1]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat1 = totalEff.hat, pers.hat1 = pers.hat)

bolb <- bol2Data[[2]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat2 = totalEff.hat, pers.hat2 = pers.hat)

bolc <- bol2Data[[3]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat3 = totalEff.hat, pers.hat3 = pers.hat)

bold <- bol2Data[[4]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat4 = totalEff.hat, pers.hat4 = pers.hat)

bole <- bol2Data[[5]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat5 = totalEff.hat, pers.hat5 = pers.hat)

# Combine all predictions
bol2_gbms <- bola %>%
  full_join(bolb) %>%
  full_join(bolc) %>%
  full_join(bold) %>%
  full_join(bole)

# Merge the predictions with the original dataset
bol2 <- merge(bol2, bol2_gbms, by = "M", all = T)

# Create the pork variable for Bolivia2
bol2$pork <- ifelse(bol2$sedepe10 > 0 | # sedepe11, 19, 36 missing
  bol2$sedepe18 > 0 |
  bol2$sedepe28, 1, 0)

bol2$pork <- replace_na(bol2$pork, 0)

# with (bol2, rowSums (cbind (pork, pork_committee)))

##################
#### Portugal ####
##################

por <- read.csv("data/ch12/csv/Portugal.csv")

# Extract all unique district magnitudes used in list elections
M <- c(unique(por$magnitude_list))

# Create the M variable (district magnitude), set to zero for non-elected candidates
por$M <- por$elected_list * por$magnitude_list

# Keep only observations with valid (non-missing, positive) magnitude
por <- por[por$M >= 1 & !is.na(por$M), ]

# Create por_real_cases Dataset

por_real_cases <- rbind(data.frame(data.frame(tballot = 0, ballot_type = "closed", M = M, formula = "dhondt", new.nvotes = "One", pool_level = "party", threshold = 0)))

# Format variables as factors with levels matching the trained GBM models
por_real_cases$new.nvotes <- factor(por_real_cases$new.nvotes, levels = optimalGBMIntra[[1]]$var.levels[[3]])
por_real_cases$pool_level <- factor(por_real_cases$pool_level, levels = optimalGBMIntra[[1]]$var.levels[[2]])
por_real_cases$ballot_type <- factor(por_real_cases$ballot_type, levels = optimalGBMIntra[[1]]$var.levels[[1]])
por_real_cases$formula <- factor(por_real_cases$formula, levels = optimalGBMIntra[[1]]$var.levels[[6]])

# Sort cases by ballot type and district magnitude
por_real_cases <- por_real_cases[order(por_real_cases$ballot_type, por_real_cases$M), ]


# Generate predictions from GBM models (both inter- and intraparty) for 5 imputed datasets
porData <- list()
for (i in 1:5) {
  # Interparty
  por_real_cases$totalEff.hat <- predict(optimalGBMInter[[i]], por_real_cases, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  por_real_cases$pers.hat <- predict(optimalGBMIntra[[i]], por_real_cases, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  porData[[i]] <- por_real_cases
}

# Extract and rename predictions for each imputation
por1 <- porData[[1]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat1 = totalEff.hat, pers.hat1 = pers.hat)

por2 <- porData[[2]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat2 = totalEff.hat, pers.hat2 = pers.hat)

por3 <- porData[[3]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat3 = totalEff.hat, pers.hat3 = pers.hat)

por4 <- porData[[4]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat4 = totalEff.hat, pers.hat4 = pers.hat)

por5 <- porData[[5]] %>%
  dplyr::select(formula, M, tballot, totalEff.hat, pers.hat, ballot_type) %>%
  rename(totalEff.hat5 = totalEff.hat, pers.hat5 = pers.hat)

# Combine all five sets of predictions into a single dataset
por_gbms <- por1 %>%
  full_join(por2) %>%
  full_join(por3) %>%
  full_join(por4) %>%
  full_join(por5)

# Merge the predictions with the original Portugal dataset by district magnitude (M)
por <- merge(por, por_gbms, by = "M", all = T)

# Create the 'pork' variable: 1 if the MP is assigned to any distributive (pork-barrel) committee
por$pork <- ifelse(por$sedepe10 > 0 | # sedepe19, 36 missing
  por$sedepe11 > 0 |
  por$sedepe18 > 0 |
  por$sedepe28 > 0, 1, 0)


por$pork <- replace_na(por$pork, 0)

# with (por, rowSums (cbind (pork, porkcommittee)))

# #########################################
# #### Competitive environment clusters ####
# ##########################################
#
# # We cluster standard errors at the "competitive environment" level
# # For our purposes, the "competitive environment" of MPs is the AP-country-year group in which they find themselves
#
# ger$comp.env <-  factor (paste (ger$pers.hat1, ger$country, ger$election_year, sep="-"))
# jpn$comp.env <-  factor (paste (jpn$pers.hat1, jpn$country, jpn$election_year, sep="-"))
# jpn1$comp.env <-  factor (paste (jpn1$pers.hat1, jpn1$country, jpn1$election_year, sep="-"))
# nzl$comp.env <-  factor (paste (nzl$pers.hat1, nzl$country, nzl$election_year, sep="-"))
# bol$comp.env <-  factor (paste (bol$pers.hat1, bol$country, bol$electionyear, sep="-"))
# por$comp.env <-  factor (paste (por$pers.hat1, por$country, por$electionyear, sep="-"))

##########################
#### Merging datasets ####
##########################

# For each country, keep only relevant variables, rename where necessary, and drop rows with missing country info

ger <- ger %>%
  dplyr::select(
    country, election_year, pork, female, incumbent, party_name,
    minister, localassembly,
    totalEff.hat1, totalEff.hat2, totalEff.hat3,
    totalEff.hat4, totalEff.hat5, pers.hat1,
    pers.hat2, pers.hat3, pers.hat4, pers.hat5, ballot_type, M
  ) %>%
  rename(year = election_year) %>%
  dplyr::filter(is.na(country) == F)

jpn <- jpn %>%
  dplyr::select(
    country, election_year, pork, female, incumbent, party_name,
    minister, localassembly,
    totalEff.hat1, totalEff.hat2, totalEff.hat3,
    totalEff.hat4, totalEff.hat5, pers.hat1,
    pers.hat2, pers.hat3, pers.hat4, pers.hat5, ballot_type, M
  ) %>%
  rename(year = election_year) %>%
  dplyr::filter(is.na(country) == F)

jpn1 <- jpn1 %>%
  dplyr::select(
    country, election_year, pork, female, incumbent, party_name,
    minister, localassembly,
    totalEff.hat1, totalEff.hat2, totalEff.hat3,
    totalEff.hat4, totalEff.hat5, pers.hat1,
    pers.hat2, pers.hat3, pers.hat4, pers.hat5, ballot_type, M
  ) %>%
  rename(year = election_year) %>%
  dplyr::filter(is.na(country) == F)

nzl <- nzl %>%
  dplyr::select(
    country, election_year, pork, female, incumbent, party_name,
    minister, localassembly,
    totalEff.hat1, totalEff.hat2, totalEff.hat3,
    totalEff.hat4, totalEff.hat5, pers.hat1,
    pers.hat2, pers.hat3, pers.hat4, pers.hat5, ballot_type, M
  ) %>%
  rename(year = election_year) %>%
  dplyr::filter(is.na(country) == F)

nzl1 <- nzl1 %>%
  dplyr::select(
    country, election_year, pork, female, incumbent, party_name,
    minister, localassembly,
    totalEff.hat1, totalEff.hat2, totalEff.hat3,
    totalEff.hat4, totalEff.hat5, pers.hat1,
    pers.hat2, pers.hat3, pers.hat4, pers.hat5, ballot_type, M
  ) %>%
  rename(year = election_year) %>%
  dplyr::filter(is.na(country) == F)

bol <- bol %>%
  dplyr::select(
    country, electionyear, pork, female, incumbent, party, minister, localassembly,
    totalEff.hat1, totalEff.hat2, totalEff.hat3,
    totalEff.hat4, totalEff.hat5, pers.hat1,
    pers.hat2, pers.hat3, pers.hat4, pers.hat5, ballot_type, M
  ) %>%
  rename(party_name = party, year = electionyear) %>%
  dplyr::filter(is.na(country) == F)

bol1 <- bol1 %>%
  dplyr::select(
    country, electionyear, pork, female, incumbent, party, minister, localassembly,
    totalEff.hat1, totalEff.hat2, totalEff.hat3,
    totalEff.hat4, totalEff.hat5, pers.hat1,
    pers.hat2, pers.hat3, pers.hat4, pers.hat5, ballot_type, M
  ) %>%
  rename(party_name = party, year = electionyear) %>%
  dplyr::filter(is.na(country) == F)

bol2 <- bol2 %>%
  dplyr::select(
    country, electionyear, pork, female, incumbent, party, minister, localassembly,
    totalEff.hat1, totalEff.hat2, totalEff.hat3,
    totalEff.hat4, totalEff.hat5, pers.hat1,
    pers.hat2, pers.hat3, pers.hat4, pers.hat5, ballot_type, M
  ) %>%
  rename(party_name = party, year = electionyear) %>%
  dplyr::filter(is.na(country) == F)

por <- por %>%
  dplyr::select(
    country, electionyear, pork, female, incumbent_election, party_name, minister, localassembly,
    totalEff.hat1, totalEff.hat2, totalEff.hat3,
    totalEff.hat4, totalEff.hat5, pers.hat1,
    pers.hat2, pers.hat3, pers.hat4, pers.hat5, ballot_type, M
  ) %>%
  rename(incumbent = incumbent_election, year = electionyear) %>%
  dplyr::filter(is.na(country) == F)

# Merge all datasets into a single object
dta <- rbind(ger, jpn, jpn1, bol, bol1, bol2, por, nzl) # , nzl1

# Handle missing data in key binary variables by setting NA to 0
dta$minister <- replace_na(dta$minister, 0)
dta$localassembly <- replace_na(dta$localassembly, 0)

# Ensure incumbent is binary and not NA
dta$incumbent <- as.numeric(dta$incumbent)
dta$incumbent <- replace_na(dta$incumbent, 0)
dta$incumbent <- as.numeric(ifelse(dta$incumbent > 0, 1, 0))

# Create a count column (1 per observation) to calculate group sizes
dta$count <- 1
dta$year <- as.factor(dta$year)
dta <- dta %>%
  group_by(country, year) %>%
  mutate(total = sum(count))

# Compute total number of MPs per country-year
dta <- dta %>%
  group_by(country, year, party_name) %>%
  mutate(total_party = sum(count))

# Calculate the proportion of seats held by each party
# 'prop' is an individual-level rendition of a legislator's
dta$prop <- dta$total_party / dta$total

# Compute number of MPs per party assigned to a pork-barrel committee
dta <- dta %>%
  group_by(country, year, party_name) %>%
  mutate(total_pork = sum(pork))

# Clean the environment: keep only the merged dataset and previously loaded country-level data
rm(list = setdiff(ls(), c("dta", "dta_merged")))

# Create a new variable that bins the predicted personalism score (pers.hat) into quartiles
dta_merged$pers.hat.quart <- cut(
  dta_merged$pers.hat,
  c(
    0,
    quantile(dta_merged$pers.hat,
      prob = c(0.25, 0.5, 0.75)
    ),
    max(dta_merged$pers.hat)
  )
)

# Set plotting area for the histogram to occupy the left 40% of the horizontal space
par(fig = c(0, 0.4, 0, 1))

# Create histogram of the number of standing committees (numscl)
hist(dta_merged$numscl,
  breaks = 120, axes = F,
  xlab = "Number of committees",
  xlim = c(0, 100),
  main = ""
)
# Add y-axis ticks
axis(2)
# Add custom x-axis ticks at regular intervals
axis(1, at = c(0, 20, 40, 80, 100), labels = c(0, 20, 40, 80, 100))

# Annotate plot with the mean number of committees
text(x = 65, y = 8, labels = paste0("mean=", round(mean(dta_merged$numscl), 2)), pos = 4, cex = 1.1)

# Annotate plot with the variance of number of committees
text(x = 65, y = 7.7, labels = paste0("var=", round(var(dta_merged$numscl), 2)), pos = 4, cex = 1.1)

# Define new plotting area for the boxplot, occupying the remaining 60% of the space
par(fig = c(.41, 1, 0, 1), new = T)

# Create boxplot of number of committees by personalism quartiles
with(
  dta_merged,
  boxplot(numscl ~ pers.hat.quart,
    xlab = "Average personalism (quartiles)",
    ylab = "Number of committees"
  )
)

# Drop extra AP and TDE columns and rename the first set as "AP" and "TDE" for model 1
smallData1 <- dta_merged %>%
  dplyr::select(-c(pers.hat2, pers.hat3, pers.hat4, pers.hat5, totalEff.hat2, totalEff.hat3, totalEff.hat4, totalEff.hat5)) %>%
  rename(AP = pers.hat1, TDE = totalEff.hat1)

# Repeat for the four other imputed datasets, rotating which AP/TDE is kept

# retains pers.hat2 and totalEff.hat2
smallData2 <- dta_merged %>%
  dplyr::select(-c(pers.hat1, pers.hat3, pers.hat4, pers.hat5, totalEff.hat1, totalEff.hat3, totalEff.hat4, totalEff.hat5)) %>%
  rename(AP = pers.hat2, TDE = totalEff.hat2)

# retains pers.hat3 and totalEff.hat3
smallData3 <- dta_merged %>%
  dplyr::select(-c(pers.hat1, pers.hat2, pers.hat4, pers.hat5, totalEff.hat1, totalEff.hat2, totalEff.hat4, totalEff.hat5)) %>%
  rename(AP = pers.hat3, TDE = totalEff.hat3)

# retains pers.hat4 and totalEff.hat4
smallData4 <- dta_merged %>%
  dplyr::select(-c(pers.hat1, pers.hat2, pers.hat3, pers.hat5, totalEff.hat1, totalEff.hat2, totalEff.hat3, totalEff.hat5)) %>%
  rename(AP = pers.hat4, TDE = totalEff.hat4)

# retains pers.hat5 and totalEff.hat5
smallData5 <- dta_merged %>%
  dplyr::select(-c(pers.hat1, pers.hat2, pers.hat3, pers.hat4, totalEff.hat1, totalEff.hat2, totalEff.hat3, totalEff.hat4)) %>%
  rename(AP = pers.hat5, TDE = totalEff.hat5)

# Combine all datasets into a list
allData <- list(smallData1, smallData2, smallData3, smallData4, smallData5)

# Remove the individual objects to clean memory
rm(smallData1, smallData2, smallData3, smallData4, smallData5)


## First model
#  Run Poisson regressions predicting the number of committees from AP and TDE (only)


lm1 <- glm(numscl ~ AP + TDE,
  family = "poisson", data = allData[[1]]
)
lm2 <- glm(numscl ~ AP + TDE,
  family = "poisson", data = allData[[2]]
)
lm3 <- glm(numscl ~ AP + TDE,
  family = "poisson", data = allData[[3]]
)
lm4 <- glm(numscl ~ AP + TDE,
  family = "poisson", data = allData[[4]]
)
lm5 <- glm(numscl ~ AP + TDE,
  family = "poisson", data = allData[[5]]
)


# Combine results using Rubin’s rules (multiple imputation)
list0 <- list(lm1, lm2, lm3, lm4, lm5)
out2 <- summary(mice::pool(list0))

# Replace them in a lm file, we will use a copy of lm3
lm_c1 <- lm3

# Get Betas and SE
se1 <- out2$std.error
for (i in 1:nrow(out2)) {
  lm_c1[["coefficients"]][[i]] <- out2$estimate[i]
}


## Second model
#  Estimate Poisson models adding institutional controls

lm1 <- glm(
  numscl ~ AP + TDE +
    v2psprlnks_osp + numberlower + executive +
    bicameral + restrictionappt + logg,
  family = "poisson", data = allData[[1]]
)
lm2 <- glm(
  numscl ~ AP + TDE +
    v2psprlnks_osp + numberlower + executive +
    bicameral + restrictionappt + logg,
  family = "poisson", data = allData[[2]]
)
lm3 <- glm(
  numscl ~ AP + TDE +
    v2psprlnks_osp + numberlower + executive +
    bicameral + restrictionappt + logg,
  family = "poisson", data = allData[[3]]
)
lm4 <- glm(
  numscl ~ AP + TDE +
    v2psprlnks_osp + numberlower + executive +
    bicameral + restrictionappt + logg,
  family = "poisson", data = allData[[4]]
)
lm5 <- glm(
  numscl ~ AP + TDE +
    v2psprlnks_osp + numberlower + executive +
    bicameral + restrictionappt + logg,
  family = "poisson", data = allData[[5]]
)


# Pool and extract results as above
list0 <- list(lm1, lm2, lm3, lm4, lm5)
out2 <- summary(mice::pool(list0))

# Replace them in a lm file, we will use a copy of lm3
lm_c2 <- lm3

# Get Betas and SE
se2 <- out2$std.error
for (i in 1:nrow(out2)) {
  lm_c2[["coefficients"]][[i]] <- out2$estimate[i]
}

### Negative binomial
# Run Negative Binomial regression to address overdispersion
# Use log(number of legislators) as an offset
#
# IMPORTANT: the code below is reproduced EXACTLY as published, including a
# mistake. `offset(log(numberlower))` is passed after the formula and a comma,
# as a separate argument. That does not create an offset.
#
# R matches named arguments first and then fills the remaining formals in
# order. Since `data` is supplied by name, the expression lands in glm.nb's
# third formal, `weights`. And offset() outside a formula is simply the
# identity function -- it returns its argument unchanged. So these models are
# WEIGHTED by log(number of legislators) rather than offset by it.
#
# It is kept as-is here so that Table 12.1 reproduces the printed table. The
# section immediately after the table sets out what the intended model gives.

lm1 <- glm.nb(
  numscl ~ AP + TDE +
    v2psprlnks_osp + executive +
    bicameral + restrictionappt + logg,
  offset(log(numberlower)),
  data = allData[[1]]
)
lm2 <- glm.nb(
  numscl ~ AP + TDE +
    v2psprlnks_osp + executive +
    bicameral + restrictionappt + logg,
  offset(log(numberlower)),
  data = allData[[2]]
)
lm3 <- glm.nb(
  numscl ~ AP + TDE +
    v2psprlnks_osp + executive +
    bicameral + restrictionappt + logg,
  offset(log(numberlower)),
  data = allData[[3]]
)
lm4 <- glm.nb(
  numscl ~ AP + TDE +
    v2psprlnks_osp + executive +
    bicameral + restrictionappt + logg,
  offset(log(numberlower)),
  data = allData[[4]]
)
lm5 <- glm.nb(
  numscl ~ AP + TDE +
    v2psprlnks_osp + executive +
    bicameral + restrictionappt + logg,
  offset(log(numberlower)),
  data = allData[[5]]
)


# Pool and replace coefficients as before
list0 <- list(lm1, lm2, lm3, lm4, lm5)
out2 <- summary(mice::pool(list0))

# Replace them in a lm file, we will use a copy of lm3
lm_c3 <- lm3

# Get Betas and SE
se3 <- out2$std.error
for (i in 1:nrow(out2)) {
  lm_c3[["coefficients"]][[i]] <- out2$estimate[i]
}


# Generate text output of all three models with stargazer
stargazer(lm_c1, lm_c2, lm_c3,
  model.names = F,
  se = list(se1, se2, se3),
  dep.var.labels = "Number of committees",
  covariate.labels = c(
    "AP", "TDE",
    "Programmatic parties",
    "Number of legislators",
    "Parliamentary regime",
    "Bicameral legislature",
    "Constitutional constraints",
    "GDP per capita (log)"
  ), type = "text"
)

# Latex table
# stargazer(lm_c1, lm_c2, lm_c3,
#                      caption = "Predictors of the number of committees in 70 legislatures observed around 2015. Models 1 and 2 are Poisson regression models; Model 3 is a negative binomial model that uses number of legislators (log scale) as an offset.",
#                      label = "T:committee_size",
#                      model.names = F,
#                      se = list(se1, se2, se3),
#                      dep.var.labels = "Number of committees",
#                      covariate.labels = c("AP", "TDE",
#                                        "Programmatic parties",
#                                        "Number of legislators",
#                                        "Parliamentary regime",
#                                        "Bicameral legislature",
#                                        "Constitutional constraints",
#                                        "GDP per capita (log)"),
#                      type="latex")

# ---------------------------------------------------------------------------
# Three versions of the Model 3 negative binomial, fitted to the same five
# datasets and pooled the same way. They differ only in how chamber size is
# handled:
#
#   (1) AS PUBLISHED   offset() outside the formula -> read as WEIGHTS.
#                      Chambers weighted by log(size); outcome stays a count.
#   (2) CORRECTED      offset() inside the formula -> a genuine offset.
#                      Outcome becomes a rate: committees per legislator.
#   (3) NO ADJUSTMENT  chamber size left out entirely; outcome is a raw count.
#
# Column (1) reproduces column 3 of Table 12.1 exactly.
# ---------------------------------------------------------------------------

# Fit one specification across the five datasets and pool with Rubin's rules,
# writing the pooled coefficients into a fitted object so stargazer can display
# them -- the same device used for Table 12.1.
pool_nb <- function(fitfun) {
  fits <- lapply(allData, fitfun)
  out <- summary(mice::pool(fits))
  obj <- fits[[3]]
  for (i in 1:nrow(out)) obj[["coefficients"]][[i]] <- out$estimate[i]
  list(model = obj, se = out$std.error)
}

# (1) As published: offset() sits outside the formula, so it becomes `weights`
nb_published <- pool_nb(function(d) {
  glm.nb(
    numscl ~ AP + TDE + v2psprlnks_osp + executive +
      bicameral + restrictionappt + logg,
    offset(log(numberlower)),
    data = d
  )
})

# (2) Corrected: offset() inside the formula, so it is a genuine offset
nb_offset <- pool_nb(function(d) {
  glm.nb(numscl ~ AP + TDE + v2psprlnks_osp + executive +
    bicameral + restrictionappt + logg +
    offset(log(numberlower)), data = d)
})

# (3) No adjustment for chamber size at all
nb_plain <- pool_nb(function(d) {
  glm.nb(numscl ~ AP + TDE + v2psprlnks_osp + executive +
    bicameral + restrictionappt + logg, data = d)
})

stargazer(nb_published$model, nb_offset$model, nb_plain$model,
  model.names = FALSE,
  se = list(nb_published$se, nb_offset$se, nb_plain$se),
  dep.var.labels = "Number of committees",
  column.labels = c("As published", "Corrected offset", "No adjustment"),
  covariate.labels = c(
    "AP", "TDE",
    "Programmatic parties",
    "Parliamentary regime",
    "Bicameral legislature",
    "Constitutional constraints",
    "GDP per capita (log)"
  ),
  type = "text"
)

# Panel (a): Plot predicted number of committees across values of AP from Poisson model (lm_c1)

gg1 <- plot_model(lm_c1,
  terms = "AP", type = "pred", show.data = F,
  colors = "gs"
) +
  geom_line(colour = "white", linewidth = 2) +
  theme_classic(base_size = 20) +
  #  scale_color_manual(color = "grey")+
  labs(
    y = "Expected number of committees",
    x = "Average Personalism",
    title = ""
  )

# ggsave ("expected_committees_4_poisson.pdf", plot=gg1, h=6, w=9)


# Panel (b): Same as above, but using negative binomial model with offset (lm_c3)

gg2 <- plot_model(lm_c3,
  terms = "AP", type = "pred", show.data = F,
  colors = "bw"
) +
  geom_line(colour = "white", linewidth = 2) +
  theme_classic(base_size = 20) +
  #  scale_color_manual(color = "grey")+
  labs(
    y = "Expected number of committees",
    x = "Average Personalism",
    title = ""
  )

# ggsave ("expected_committees_6_negbin.pdf", plot=gg2, h=6, w=9)

# Combine both plots in a single row layout
grid.arrange(gg1, gg2, ncol = 2)

# ---------------------------------------------------------------------------
# Rebuild the chamber-level data from scratch, using the LATER scores file.
#
# Everything is rebuilt rather than reused because the objects from Part I are
# no longer in memory -- an rm() further up cleared all but `dta` and
# `dta_merged`, and `dta` has since been reassigned to the country-level data
# used by Part II. New names (suffix B) are used throughout so that nothing
# Part II relies on is disturbed.
# ---------------------------------------------------------------------------

# The legislative-institutions data, read again under a fresh name.
vdatB <- read_dta("data/ch12/dta/VanDuskyAllenTouchtonPRQ.dta")

dta_smallB <- vdatB %>%
  dplyr::select(
    year, country_name, numscl, v2psprlnks_osp, numberlower,
    eleclower, executive, effn, bicameral, restrictionappt, logg
  ) %>%
  rename(country = country_name)

dta_smallB$country <- dplyr::recode(dta_smallB$country,
  "Korea South" = "South Korea",
  "Slovak Republic" = "Slovakia",
  "United States" = "United States of America"
)

# Load the August 2024 scores into their OWN environment, so the object named
# data2export that Part I created is not overwritten.
envB <- new.env()
load("data/shared/RealSystems_Scores_GBM_Aug_2024.RData", envir = envB)

# Same 14 columns as Part I. This matters here: the August 2024 file carries ten
# extra columns, three of them with missing values, and drop_na() below would
# otherwise discard observations on the strength of variables never used.
scoreColsB <- c(
  "country", "year",
  paste0("totalEff.hat", 1:5),
  paste0("pers.hat", 1:5),
  "totalEff.hat", "pers.hat"
)

# Merge, filter and clean exactly as in Part I.
dta_mergedB <- merge(envB$data2export[, scoreColsB], dta_smallB,
  by = "country", all.y = TRUE
)
dta_mergedB <- dta_mergedB[dta_mergedB$year.x >= 2005 & dta_mergedB$year.x <= 2015, ]
dta_mergedB <- dta_mergedB %>%
  drop_na() %>%
  arrange(country, year.x)

# Keep only the election immediately prior to 2015 for each country.
dta_mergedB$proper_elections <- 0
for (i in 1:(nrow(dta_mergedB) - 1)) {
  if (dta_mergedB$country[i] == dta_mergedB$country[i + 1]) {
    dta_mergedB$proper_elections[i] <- 1
  }
}
dta_mergedB <- dta_mergedB %>%
  dplyr::filter(proper_elections != 1) %>%
  dplyr::select(-proper_elections)

# Five datasets, rotating which imputed AP/TDE pair is treated as the predictor.
# Written as a loop rather than five copies, but the result is the same.
allDataB <- lapply(1:5, function(i) {
  keepAP <- paste0("pers.hat", i)
  keepTDE <- paste0("totalEff.hat", i)
  drop <- c(
    setdiff(paste0("pers.hat", 1:5), keepAP),
    setdiff(paste0("totalEff.hat", 1:5), keepTDE)
  )
  x <- dta_mergedB[, !(names(dta_mergedB) %in% drop)]
  names(x)[names(x) == keepAP] <- "AP"
  names(x)[names(x) == keepTDE] <- "TDE"
  x
})

# ---------------------------------------------------------------------------
# The same three models as Table 12.1, with the same specifications.
# Each is fitted to all five datasets and pooled with Rubin's rules, then the
# pooled coefficients are written into one fitted object so stargazer can
# display them -- the same device used in Part I.
# ---------------------------------------------------------------------------

fit_and_pool <- function(fitfun) {
  fits <- lapply(allDataB, fitfun)
  out <- summary(mice::pool(fits))
  obj <- fits[[3]] # a copy to carry the pooled numbers
  for (i in 1:nrow(out)) obj[["coefficients"]][[i]] <- out$estimate[i]
  list(model = obj, se = out$std.error)
}

# Model 1: AP and TDE only
m1B <- fit_and_pool(function(d) {
  glm(numscl ~ AP + TDE, family = "poisson", data = d)
})

# Model 2: adding institutional controls
m2B <- fit_and_pool(function(d) {
  glm(numscl ~ AP + TDE + v2psprlnks_osp + numberlower + executive +
    bicameral + restrictionappt + logg, family = "poisson", data = d)
})

# Model 3: negative binomial, as in Part I
# Note this reproduces Part I's specification exactly, mistake included: the
# offset term sits outside the formula and is therefore read as weights. That
# is deliberate. The purpose of this section is to isolate the effect of
# changing the SCORES FILE, so every other feature of the model must be held
# constant. The separate question of the offset is dealt with in the correction
# section that follows Table 12.1.
m3B <- fit_and_pool(function(d) {
  glm.nb(
    numscl ~ AP + TDE + v2psprlnks_osp + executive +
      bicameral + restrictionappt + logg,
    offset(log(numberlower)),
    data = d
  )
})

cat(
  "Observations:", nrow(dta_mergedB),
  " (Part I, on the earlier scores, used 70)\n\n"
)

stargazer(m1B$model, m2B$model, m3B$model,
  model.names = FALSE,
  se = list(m1B$se, m2B$se, m3B$se),
  dep.var.labels = "Number of committees",
  covariate.labels = c(
    "AP", "TDE",
    "Programmatic parties",
    "Number of legislators",
    "Parliamentary regime",
    "Bicameral legislature",
    "Constitutional constraints",
    "GDP per capita (log)"
  ),
  type = "text"
)

# Start by constructing five datasets, one per imputed AP/TDE draw

smallData1 <- dta %>%
  dplyr::select(-c(pers.hat2, pers.hat3, pers.hat4, pers.hat5, totalEff.hat2, totalEff.hat3, totalEff.hat4, totalEff.hat5)) %>%
  rename(AP = pers.hat1, TDE = totalEff.hat1)

smallData2 <- dta %>%
  dplyr::select(-c(pers.hat1, pers.hat3, pers.hat4, pers.hat5, totalEff.hat1, totalEff.hat3, totalEff.hat4, totalEff.hat5)) %>%
  rename(AP = pers.hat2, TDE = totalEff.hat2)

smallData3 <- dta %>%
  dplyr::select(-c(pers.hat1, pers.hat2, pers.hat4, pers.hat5, totalEff.hat1, totalEff.hat2, totalEff.hat4, totalEff.hat5)) %>%
  rename(AP = pers.hat3, TDE = totalEff.hat3)

smallData4 <- dta %>%
  dplyr::select(-c(pers.hat1, pers.hat2, pers.hat3, pers.hat5, totalEff.hat1, totalEff.hat2, totalEff.hat3, totalEff.hat5)) %>%
  rename(AP = pers.hat4, TDE = totalEff.hat4)

smallData5 <- dta %>%
  dplyr::select(-c(pers.hat1, pers.hat2, pers.hat3, pers.hat4, totalEff.hat1, totalEff.hat2, totalEff.hat3, totalEff.hat4)) %>%
  rename(AP = pers.hat5, TDE = totalEff.hat5)


# Create a party-country identifier to allow nested analysis or fixed effects
smallData1$party_name_country <- paste(smallData1$party_name, smallData1$country, sep = "-")
smallData2$party_name_country <- paste(smallData2$party_name, smallData2$country, sep = "-")
smallData3$party_name_country <- paste(smallData3$party_name, smallData3$country, sep = "-")
smallData4$party_name_country <- paste(smallData4$party_name, smallData4$country, sep = "-")
smallData5$party_name_country <- paste(smallData5$party_name, smallData5$country, sep = "-")


# Combine all datasets into one list and clear memory
allData <- list(smallData1, smallData2, smallData3, smallData4, smallData5)

rm(smallData1, smallData2, smallData3, smallData4, smallData5)


# Create five boxplots displaying variation of AP per country
# g <- list()
# for (i in 1:5){
#   g[[i]] <- ggplot(allData[[i]], aes(x = country, y = AP)) +
#     geom_boxplot() +
#     ggtitle(paste("Boxplot of AP by Country, sample", i, sep=" ")) +
#     xlab("Country") +
#     ylab("AP")
# }
# grid.arrange(g[[1]],g[[2]],g[[3]],g[[4]],g[[5]], ncol=3)