How to read this file

This is a replication file. It contains every line of R code needed to reproduce the tables and figures of Chapter 6, together with an explanation of what each line does. The code is exactly the code we ran for the book; nothing has been simplified.

Chapter 6 employs gradient-boosting machine (GBM) — “machine learning” — models to sift through the enormous amount of information contained in the simulated worlds of Chapter 4, in order to identify the best set of predictors for TDE and AP scores. That knowledge is then combined with what we know about specific elections — the 2013 German federal election, say — to produce scores for the interparty and intraparty incentives generated by their component rules. We leverage a comprehensive collection of detailed data on real elections held in 140 countries between 1945 and 2015, recording each electoral component used in Chapter 4, and use them to generate TDE and AP for all the systems used in 1,528 elections.

The shape of this file. Chapter 6 is the longest replication file in the book, and it runs in three stages:

A worked example (Tables 6.1 and 6.2). Seven districts in France, Bulgaria and Spain, used to show concretely how a description of an electoral system becomes a predicted TDE and AP score, and how district-level scores can be aggregated to the country level three different ways.
The full dataset (the long code-real-cases and predictions-real-elections blocks). The same operation applied to every election in the collection, plus the bookkeeping needed to label 140 countries consistently.
The figures (6.1 through 6.5). Projections of real systems onto the I–I space.

If you are reading to understand the method, stage 1 is the part to study; stage 2 is mostly careful data handling.

If you have never used R. Click the grey Code button above any block to show or hide the code. Every line beginning with # is a comment, ignored by R and written purely for you. A few conventions used throughout:

<- is R’s assignment arrow; x <- 5 stores 5 under the name x.
library(name) loads an add-on package of extra functions.
list() holds several objects at once; mylist[[3]] retrieves the third.
for (i in 1:5) { ... } repeats the braced instructions five times.
df$col pulls the column col out of the data frame df.
%>% is the “pipe”: a %>% f() is another way of writing f(a).

Setup

This chapter loads a long list of packages. Not all are used by the code that remains in this file — several are carried over from the model-fitting stage documented in the companion simulation file. If R reports there is no package called ..., install it once with install.packages("packagename").

Why dplyr:: is spelled out in the code below. filter() and select() have namesakes in other packages — filter() in base R’s stats, select() in MASS, and this chapter also loads car. R resolves a name to whichever package was attached most recently, so a bare filter() can silently become the wrong function. Writing dplyr::filter() removes the ambiguity. Note that library(tidyverse) is no defence: loading an already-attached package does nothing, so it cannot restore dplyr’s precedence.

# Load relevant libraries

# The workhorses for this chapter are gbm (the fitted prediction models),
# ggplot2 (every figure), and the dplyr/tidyr/tidyverse family (data handling).
# The remainder support modelling and plotting steps used elsewhere in the
# project: mgcv and rpart fit alternative model types, caret manages model
# training, mixtools and ellipse draw the
# confidence ellipses in the I–I space figures, and doMC enables parallel
# processing.
library(mgcv)
library(ggplot2)
library(gridExtra)
library(rpart)
library(car)
library(gbm)
library(caret)
library(dplyr)
library(tidyr)
library(magrittr)
library(tidyverse)
require(mixtools)
library(ellipse)
library(gtools)
library(ggpubr)
library(doMC)
library(openxlsx)
library(dplyr)

Introduction

Chapter 6 builds on the theoretical and simulation foundations established in previous chapters to place real-world electoral systems within the Interparty–Intraparty (I–I) space. While the earlier chapters developed measures of electoral incentives—the Total Duvergerian Effect (TDE) capturing interparty competition, and Average Personalism (AP) capturing intraparty competition—this chapter advances the project by applying predictive models to actual electoral systems from around the world. Using gradient-boosted machine (GBM) models trained on extensive simulations, we estimate TDE and AP scores for over 1,500 elections held in 140 countries between 1945 and 2015.

This approach allows us to generate precise expectations about how the component rules of each electoral system—such as district magnitude, ballot type, vote pooling level, seat allocation formula, and electoral thresholds—combine to shape the incentives for interparty coordination and intraparty competition in real electoral contexts. By projecting actual electoral systems onto the I–I space, we identify the “ideal types” that correspond to real electoral rules and analyze their implications for political behavior. Furthermore, the chapter introduces a user-friendly interactive dashboard (the I3 dashboard) that enables researchers and practitioners to explore how variations in electoral rules influence the inter- and intraparty incentives.

This file contains calls to datasets used in Chapter 6 (Placing “Real-World” Electoral Systems in the I–I Space), as well as the code necessary to produce all graphs in the chapter.

Paths in this replication package. Every data file is addressed relative to the folder holding this document, as data/shared/... for files several chapters use, or data/ch06/... and so on for files specific to one chapter.

Nothing needs editing. knitr runs each code block with the working directory set to the document’s own folder, which is the root of this repository, so the paths resolve wherever the repository is cloned. There are no setwd() calls and no absolute paths anywhere in these files.

# resolves the short file names used below. See the note above on why this call
# is repeated in later blocks rather than stated once here.

We saved a separate .Rdata file with information on Bulgarian, French, and Spanish districts. Notice that the only factor that really varies across districts is M. The following chunk imports this spreadsheet.

# Load data containing real district-level electoral rules for selected elections in France, Bulgaria, and Spain

load("data/ch06/RData/bul_fra_spa_district_mag.RData")

The next chunk generates predictions of TDE and AP scores based on district-level GBM models.

This block completes the district-level prediction workflow described in Table 6.1 of the book, using trained models to assign predicted values of interparty (TDE) and intraparty (AP) incentives for each real district.

What a GBM prediction is doing here. A gradient boosting machine is a prediction engine built from hundreds of small decision trees, each correcting the errors of the ones before it. It was trained in the companion simulation file on the simulated worlds of Chapter 4, where the true TDE and AP are known because we generated them.

Nothing is being estimated in this file. The models are already fitted; here they are simply asked a question: given this real district’s magnitude, ballot type, formula, pooling level and threshold, what TDE and AP would a simulated system with those same rules have produced? That is what makes it possible to place a real election in the I–I space.

As elsewhere in the book, there are five fitted models rather than one, so every prediction is made five times and averaged. The spread across those five is reported as a standard deviation in Table 6.1 and is a measure of model uncertainty.

## Using the district-level GBM object ##

# Load saved GBM model objects for TDE and AP (district-level)
# load() restores objects under the names they were saved with -- here
# totalEffENP.objects (for TDE) and totalAP.objects (for AP).
load(file = "data/shared/TDE_district_objects_t485_d9.RData")
load(file = "data/shared/AP_district_objects_t285_d13.RData") # Ensure this is the most updated version

# Rename objects to clearly distinguish as district-level
# Country-level versions of the same objects are loaded further down under the
# identical names, so renaming now prevents one from overwriting the other.
totalAP.objects.district <- totalAP.objects
totalEffENP.objects.district <- totalEffENP.objects

# Remove original variables to avoid confusion
# rm() deletes the originals, so a later load() cannot silently shadow them.
rm(totalAP.objects, totalEffENP.objects)

# Initialize lists to store predicted values for each simulation replicate
# Chained assignment: both names are set to an empty list in one statement.
totalEff.hat <- pers.hat <- list()

# Loop over the five GBM models
for (i in 1:5) {
  # Retrieve i-th GBM model for AP and TDE
  # $optimalGBM is the component holding the best-performing fitted model from
  # each of the five training runs; [[i]] takes the i-th of them.
  optimalGBMIntra <- totalAP.objects.district$optimalGBM[[i]]
  optimalGBM <- totalEffENP.objects.district$optimalGBM[[i]]

  # Predict TDE score for real district cases
  # predict(model, newdata, n.trees) applies a fitted model to fresh data.
  # `realCases` holds the real-world districts; n.trees says how many of the
  # boosted trees to use -- here all of those chosen as optimal in training.
  totalEff.hat[[i]] <- predict(optimalGBM, realCases,
    n.trees = optimalGBM$n.trees
  )

  # Predict AP score for real district cases
  pers.hat[[i]] <- predict(optimalGBMIntra, realCases,
    n.trees = optimalGBMIntra$n.trees
  )
}

# Store predictions from all five GBM models in realCases dataframe
# Ten new columns are attached to the data: five TDE predictions and five AP
# predictions, one per fitted model. Writing them out one at a time is more
# verbose than a loop but leaves the column names explicit and easy to trace.
realCases$totalEff.hat1.district <- totalEff.hat[[1]]
realCases$totalEff.hat2.district <- totalEff.hat[[2]]
realCases$totalEff.hat3.district <- totalEff.hat[[3]]
realCases$totalEff.hat4.district <- totalEff.hat[[4]]
realCases$totalEff.hat5.district <- totalEff.hat[[5]]

realCases$pers.hat1.district <- pers.hat[[1]]
realCases$pers.hat2.district <- pers.hat[[2]]
realCases$pers.hat3.district <- pers.hat[[3]]
realCases$pers.hat4.district <- pers.hat[[4]]
realCases$pers.hat5.district <- pers.hat[[5]]

# Compute the average predicted TDE across the five models
# with(data, expr) evaluates expr inside the data frame, so the columns can be
# named directly instead of writing realCases$... five times over. The result is
# the pooled prediction reported throughout the chapter.
realCases$totalEff.hat.district <- with(realCases, (totalEff.hat1.district +
  totalEff.hat2.district +
  totalEff.hat3.district +
  totalEff.hat4.district +
  totalEff.hat5.district) / 5)

# Compute the average predicted AP across the five models
realCases$pers.hat.district <- with(realCases, (pers.hat1.district +
  pers.hat2.district +
  pers.hat3.district +
  pers.hat4.district +
  pers.hat5.district) / 5)

Table 6.1: TDE and AP Scores for Seven Districts in France, Bulgaria, & Spain

This table presents predicted Total Duvergerian Effect (TDE) and Average Personalism (AP) scores for seven selected districts across France (1986), Bulgaria (2014), and Spain (2015). The districts vary considerably in district magnitude (M), from as low as 2 in Alpes-de-Haute-Provence to as high as 16 in Bouches-du-Rhône and Sofia-23, providing a useful range to examine how district size influences electoral incentives.

# Define countries and corresponding district magnitudes to include in Table 6.1
# These two vectors work as a pair, read position by position: the first entry
# of each gives France with M = 2, the second France with M = 16, and so on for
# all seven districts.
country2find <- c(
  "France", "France",
  "Bulgaria", "Bulgaria", "Bulgaria",
  "Spain", "Spain"
)

distrM2find <- c(2, 16, 4, 14, 16, 6, 11)

# Initialize empty object to store table rows
# c() creates an empty object that rbind() will grow one row at a time.
table.6.1.data <- c()


# Loop over each selected district
for (i in 1:length(country2find)) {
  cn <- country2find[i] # Country name
  dm <- distrM2find[i] # District magnitude

  # Filter row in realCases dataset for that country and district magnitude
  # Square brackets index as [rows, columns]. The condition before the comma
  # keeps rows where BOTH tests pass (& means "and"); the empty slot after the
  # comma keeps all columns.
  dat <- realCases[realCases$Country == cn & realCases$M == dm, ]
  dat <- dat[1, ] # Take only the first match (in case of duplicates)

  # Extract average predicted TDE score for district
  # grep() searches the column names for a pattern and returns the positions
  # that match, so the column is found by name rather than by counting.
  tde.m <- round(dat[, grep("totalEff.hat.district", colnames(dat))], 3)

  # Calculate standard deviation across the 5 model predictions (exclude the mean column)
  # The looser pattern "totalEff.hat" matches all six TDE columns -- the five
  # individual predictions plus the average. [, -6] drops the sixth, the
  # average, so the standard deviation is taken across the five models only.
  # Including the average would understate the spread.
  tde.sd <- round(sd(dat[, grep("totalEff.hat", colnames(dat))][, -6]), 2)

  # Extract average predicted AP score for district
  ap.m <- round(dat[, grep("pers.hat.district", colnames(dat))], 3)

  # Calculate standard deviation across the 5 model predictions (exclude the mean column)
  ap.sd <- round(sd(dat[, grep("pers.hat", colnames(dat))][, -6]), 2)

  # Combine results into a row
  dt.return <- data.frame(cn, dm, tde.m, tde.sd, ap.m, ap.sd)

  # Append row to final table
  # Each pass adds one row to the bottom of the accumulating table.
  table.6.1.data <- rbind(table.6.1.data, dt.return)
}


# Print Table 6.1 data
print(table.6.1.data)

##         cn dm tde.m tde.sd  ap.m ap.sd
## 1   France  2 1.295   0.14 0.418  0.01
## 2   France 16 0.963   0.07 0.399  0.00
## 3 Bulgaria  4 0.964   0.13 0.508  0.01
## 4 Bulgaria 14 0.943   0.09 0.495  0.01
## 5 Bulgaria 16 0.945   0.07 0.496  0.01
## 6    Spain  6 1.216   0.18 0.403  0.01
## 7    Spain 11 1.205   0.13 0.399  0.00

How to read Table 6.1. Six columns: country, district magnitude, then the mean and standard deviation of predicted TDE, then the same pair for AP. The standard deviations are not sampling error in the usual sense — they measure how much the five fitted GBM models disagree with one another. Small values mean the prediction is robust to which model produced it.

Read down the magnitude column and across to TDE: the expectation is that TDE falls as districts get larger. Then check AP across the same rows. These systems share broadly similar ballot structures, so AP should move much less — which is the point of showing both columns side by side.

Reported values in Table 6.1 are the average predictions from five gradient boosting models, with standard deviations in parentheses reflecting model uncertainty. As expected, TDE scores generally decrease with increasing district magnitude, consistent with the idea that larger districts ease constraints on party competition and entry. AP scores, by contrast, exhibit relatively little variation, suggesting that intraparty incentives remain fairly stable across these districts given similar ballot structures.

This district-level analysis illustrates how specific institutional features—particularly district magnitude—translate into measurable differences in interparty and intraparty competition.

Table 6.2

Table 6.2 compares three approaches to projecting Total Duvergerian Effect (TDE) and Average Personalism (AP) scores at the country level for France (1986), Bulgaria (2014), and Spain (2015).

First, we project TDE and AP scores for Bulgaria, France, and Spain in their entirety, based on district-weighed projections. These results are the first column in Table 6.2.

The aggregation problem. A country’s electoral system is not one system but many: each district has its own magnitude and so its own incentives. Turning district scores into a single national score requires a choice, and there is no uniquely correct answer. Table 6.2 therefore reports three, and this block produces the first:

M/L weighted average (this block) — average the district scores, weighting each district by its share of the chamber’s seats. A district electing 16 of 500 seats counts sixteen times as much as one electing a single seat.
District projection (next block, column 2) — imagine a country made up entirely of average-sized districts, and predict for that.
Country-level GBM (next block, column 3) — models trained directly on whole-country simulations, bypassing districts altogether.

The three disagree, and the size of that disagreement is itself a finding.

### For country-level GBM predictions, we need to combine M ###
# We need a weighted average, where each district's M furnishes the weight
# If L is the total number of seats in the house (summing over M throughout all districts),
# a district's weight should be M over L

# Compute district-level weights based on relative magnitude M / L (total seats)
# group_by() + mutate() is the dplyr idiom for "do this separately within each
# group, but keep every row". Because the data are grouped by Country, sum(M)
# is that country's total seats rather than the total across all countries --
# so M.weight is exactly the M/L share described above.
realCases <- realCases %>%
  group_by(Country) %>%
  mutate(M.weight = M / sum(M))


# Aggregate TDE and AP using weighted average for each country
# Group by 'country' and calculate country-level 'TDE'
# summarise() differs from mutate(): it COLLAPSES each group to a single row.
# sum(M.weight * score) is the weighted average, since the weights sum to 1.
# All ten columns -- five TDE and five AP -- are aggregated in one pass, keeping
# the five model replicates separate so their spread can be reported below.
country_level <- realCases %>%
  group_by(Country) %>%
  summarise(
    TDE1 = sum(M.weight * totalEff.hat1.district),
    TDE2 = sum(M.weight * totalEff.hat2.district),
    TDE3 = sum(M.weight * totalEff.hat3.district),
    TDE4 = sum(M.weight * totalEff.hat4.district),
    TDE5 = sum(M.weight * totalEff.hat5.district),
    AP1 = sum(M.weight * pers.hat1.district),
    AP2 = sum(M.weight * pers.hat2.district),
    AP3 = sum(M.weight * pers.hat3.district),
    AP4 = sum(M.weight * pers.hat4.district),
    AP5 = sum(M.weight * pers.hat5.district)
  )


# Format results: compute mean and SD across the 5 replicate scores
table.6.2.mlweighted <- c()
for (i in 1:nrow(country_level)) {
  dat <- country_level[i, ]
  tde.m <- round(mean(as.numeric(dat[, grep("TDE", colnames(dat))])), 2)
  tde.sd <- round(sd(dat[, grep("TDE", colnames(dat))]), 2)
  ap.m <- round(mean(as.numeric(dat[, grep("AP", colnames(dat))])), 3)
  ap.sd <- round(sd(dat[, grep("AP", colnames(dat))]), 2)
  dt.return <- data.frame(dat$Country, tde.m, tde.sd, ap.m, ap.sd)
  table.6.2.mlweighted <- rbind(table.6.2.mlweighted, dt.return)
}

# Print weighted average predictions
print(table.6.2.mlweighted)

##   dat.Country tde.m tde.sd  ap.m ap.sd
## 1    Bulgaria  0.96   0.09 0.500  0.01
## 2      France  1.00   0.09 0.403  0.00
## 3       Spain  1.25   0.15 0.404  0.00

Finally, the next snippet projects TDE and AP scores based on district-level and country-level objects. These projections are the ones that we provide for general use. The projections appear in columns 2 and 3 of Table 6.2.

# Report additional district-projection and country-GBM scores
# Load precomputed predictions using the average-district method and the country-level GBM models
# ---------------------------------------------------------------------------
# This chapter previously read RealSystems_Scores_GBM.RData, a 2023-era file.
# Chapters 7 to 10 read the later Aug 2024 file, and the split was an oversight
# rather than a decision -- the newer file had in fact already been copied into
# this chapter's own Datasets/RData folder, but the code kept pointing at the
# older one in a different subfolder.
#
# Switching was verified to be costless: for the three cases used in Table 6.2
# (France 1986, Bulgaria 2014, Spain 2015) both the district-projection and the
# country-GBM figures are identical under either file, to three decimals.
#
# Both scores files live in data/shared/ because several chapters read them.
# ---------------------------------------------------------------------------
load("data/shared/RealSystems_Scores_GBM_Aug_2024.RData")

# Filter only cases corresponding to France 1986, Bulgaria 2014, and Spain 2015
Cases <- data2export %>%
  dplyr::filter(
    country %in% c("France", "Bulgaria", "Spain"),
    year %in% c(1986, 2014, 2015)
  )


# Compute average and SD of predictions from district-level GBM projections
# The two loops below are near-identical; the only difference is the sign in
# front of grep(".district", ...). Here the columns whose names CONTAIN
# ".district" are selected -- the district-projection method. In the second loop
# the minus sign drops those same columns, leaving the country-GBM ones.
# As in Table 6.1, [, -6] excludes the averaged column so the standard
# deviation reflects the five models only.
table.6.2.district.projection <- c()
for (i in 1:nrow(Cases)) {
  dat <- Cases[i, grep(".district", colnames(Cases))]
  tde.d.m <- round(mean(as.numeric(dat[, grep("totalEff.hat", colnames(dat))][, -6])), 3)
  tde.d.sd <- round(sd(dat[, grep("totalEff.hat", colnames(dat))][, -6]), 3)
  ap.d.m <- round(mean(as.numeric(dat[, grep("pers.hat", colnames(dat))][, -6])), 3)
  ap.d.sd <- round(sd(dat[, grep("pers.hat", colnames(dat))][, -6]), 2)
  dt.return <- data.frame(Cases[i, ]$country, tde.d.m, tde.d.sd, ap.d.m, ap.d.sd)
  table.6.2.district.projection <- rbind(table.6.2.district.projection, dt.return)
}


# Compute predictions from country-level GBM models (direct estimates)
table.6.2.country.projection <- c()
for (i in 1:nrow(Cases)) {
  dat <- Cases[i, -grep(".district", colnames(Cases))]
  tde.c.m <- round(mean(as.numeric(dat[, grep("totalEff.hat", colnames(dat))][, -6])), 3)
  tde.c.sd <- round(sd(dat[, grep("totalEff.hat", colnames(dat))][, -6]), 3)
  ap.c.m <- round(mean(as.numeric(dat[, grep("pers.hat", colnames(dat))][, -6])), 3)
  ap.c.sd <- round(sd(dat[, grep("pers.hat", colnames(dat))][, -6]), 2)
  dt.return <- data.frame(Cases[i, ]$country, tde.c.m, tde.c.sd, ap.c.m, ap.c.sd)
  table.6.2.country.projection <- rbind(table.6.2.country.projection, dt.return)
}


# Print district-level projection (Column 2 of Table 6.2)
print(table.6.2.district.projection)

##   Cases.i....country tde.d.m tde.d.sd ap.d.m ap.d.sd
## 1           Bulgaria   0.966    0.104  0.501    0.01
## 2             France   0.968    0.114  0.401    0.01
## 3              Spain   1.219    0.154  0.403    0.01
## 4              Spain   1.219    0.154  0.403    0.01

# Print country-level GBM projection (Column 3 of Table 6.2)
print(table.6.2.country.projection)

##   Cases.i....country tde.c.m tde.c.sd ap.c.m ap.c.sd
## 1           Bulgaria   0.503    0.071  0.476    0.01
## 2             France   0.618    0.042  0.428    0.01
## 3              Spain   0.670    0.090  0.422    0.01
## 4              Spain   0.670    0.090  0.422    0.01

How to read Table 6.2. The table is printed here as three separate blocks — the M/L weighted average from the previous chunk, then these two — which the book assembles into a single three-column table. Read across the columns for one country at a time and ask how far apart the three methods land.

Expect the M/L weighted figures to give the highest TDE, because small districts are highly constraining and the weighting still gives them a voice. Expect the country-level GBM to give the lowest, because it models the chamber as a whole, where fragmentation is easier than in any single district. AP should be more stable across the three, since it depends mostly on ballot structure rather than on how districts are aggregated. The gap between columns is a reminder that “the electoral system of a country” is a summary, and summaries involve choices.

Table 6.2 compares three different methods for projecting Total Duvergerian Effect (TDE) and Average Personalism (AP) scores at the country level for France (1986), Bulgaria (2014), and Spain (2015). The first column presents a weighted average of district-level predictions, where each district’s influence is proportional to its legislative weight (M/L). The second column uses a district-level GBM model to simulate a hypothetical country made up entirely of average-sized districts. The third column reports results from a country-based GBM model, trained on national-level simulations. TDE scores are generally highest under the M/L-weighted approach, which captures the influence of small, highly constraining districts. In contrast, the country-based GBM projections yield much lower TDE values, reflecting the reality that national-level systems often permit greater fragmentation than any individual district would suggest. AP scores show smaller variation across methods, but country-based projections are slightly higher on average, indicating that national-level conditions may foster somewhat more personalistic electoral behavior than district-level estimates alone would imply.

Data for figures

Everything up to this point was a worked example on seven districts. From here the same operation is applied to the full collection of 1,528 elections.

This block wipes the workspace. rm(list = ls()) deletes every object created so far, including realCases and the tables just built. That is deliberate — the second half of the chapter rebuilds realCases from the full spreadsheet, and clearing first guarantees no stale object survives to be silently reused. The practical consequence is that the file must be run in order, from the top; running a later block on its own will fail because the objects it expects no longer exist.

# We start by removing all objects in *R*, then re-importing the GBM objects.
# ls() lists everything in memory and rm() removes it. See the note above.
rm(list = ls())

# Load district-level GBM models for TDE and AP
load(file = "data/shared/TDE_district_objects_t485_d9.RData")
load(file = "data/shared/AP_district_objects_t285_d13.RData")

# Rename for clarity
totalAP.objects.district <- totalAP.objects
totalEffENP.objects.district <- totalEffENP.objects

# Remove original object names to avoid ambiguity
rm(totalAP.objects, totalEffENP.objects)

# Load country-level GBM models for TDE and AP
# "nofamily" records a modelling choice: these models do not use electoral
# system family as a predictor, so a system is characterised only by its
# component rules. They are left under their original names, which is why the
# district-level pair had to be renamed above.
load(file = "data/shared/TDE_country_nofamily_objects.RData")
load(file = "data/shared/AP_country_nofamily_objects.RData")

We import a dataset with all electoral systems. Just as we did with Bulgaria, France, and Spain, we need to work on this file to make all variables comformable with the way in which our GBM objects are coded. We take this preliminary step to then predict TDE and AP scores for all electoral systems. These are the scores that we will then aggregate and discuss in the various figures of Chapter 6.

Why this block is so long. The GBM models were trained on simulated systems, where every rule was set by us and took a tidy, known value. Real election data are recorded by human coders across 140 countries and seven decades, so the same rule appears under many spellings — dHondt, D'Hondt, Hagenbach-Bischoff quota with largest remainders.

A fitted model will not accept a category it has never seen. So before any prediction can be made, every real-world label has to be mapped onto the vocabulary the models were trained with. That translation is essentially all this block does: load the spreadsheet, standardise missing values, collapse dozens of recorded labels into the handful of simulated categories, and drop the few systems that were never simulated.

It is bookkeeping, not analysis — but it is where a replication most often goes wrong, which is why it is written out explicitly rather than hidden in a helper function.

# Load the dataset of real electoral systems
# read.xlsx() reads one sheet of an Excel workbook into a data frame.
realCases <- read.xlsx("data/ch06/xlsx/Dataset_05_27_21.xlsx", sheet = "data")
# table() counts how often each value occurs; names() then lists the distinct
# values. Printed here purely to inspect what spellings the file contains --
# the result is displayed, not stored.
names(table(realCases$bg_tier1_formula)) # Check distinct formula names in tier 1 for inspection

##  [1] "-88"                                                     
##  [2] "-99"                                                     
##  [3] "Absolute Majority"                                       
##  [4] "Alternative Vote"                                        
##  [5] "Block Vote"                                              
##  [6] "dHondt"                                                  
##  [7] "Droop quota"                                             
##  [8] "Droop quota with largest remainders"                     
##  [9] "Fortified PR"                                            
## [10] "Hagenbach-Bischoff quota"                                
## [11] "Hagenbach-Bischoff quota with highest average remainders"
## [12] "Hagenbach-Bischoff quota with largest remainders"        
## [13] "Hare quota"                                              
## [14] "Hare quota with highest average remainders"              
## [15] "Hare quota with largest remainders"                      
## [16] "Imperiali quota"                                         
## [17] "Imperialli quota with largest remainders"                
## [18] "Limited Nomination"                                      
## [19] "Limited Vote"                                            
## [20] "Modified Sainte-Lague"                                   
## [21] "Party Block Vote"                                        
## [22] "Reinforced Imperiali quota"                              
## [23] "Sainte-Lague"                                            
## [24] "Single Nontransferable Vote"                             
## [25] "Single Transferable Vote"                                
## [26] "Single-Member-District-Plurality"                        
## [27] "Two-Round Absolute Majority"                             
## [28] "Two-Round Majority Runoff"                               
## [29] "Two-Round Majority-Plurality"                            
## [30] "Two-Round Party Block Vote"                              
## [31] "Two-Round Qualified Majority"

## Clean and re-code variables to match model structure
# Force number of votes to "One" for mixed-member systems
# %in% tests membership in a list. Read the line as: for rows whose electoral
# rule is one of these two, set votesN to "One".
realCases$votesN[realCases$bg_elecrule %in% c("Mixed Dependent", "Mixed Independent")] <- "One"

# Convert special missing codes to NA
# Datasets often encode "missing" with sentinel numbers rather than a blank.
# Left unconverted, -99 would be treated as a real quantity and silently
# corrupt every calculation downstream. NA is R's proper missing-value marker.
realCases[realCases == -99] <- NA
realCases[realCases == -88] <- NA

# Remove non-simulated case: Mali's 2-round Party Block Vote
# The models cannot predict for a system that was never simulated, so the case
# is dropped. != means "not equal to".
realCases <- realCases[realCases$bg_tier1_formula != "Two-Round Party Block Vote", ]


## Generate numeric threshold variable
# An alternative way to code Threshold
# Create numeric thresholds for tier 1 and 2
# Many systems have two tiers, each potentially with its own legal threshold.
# The simulations used a single threshold, so the two must be reduced to one.
t1 <- as.numeric(realCases$thresholdPerc_tier1)
t2 <- as.numeric(realCases$thresholdPerc_tier2)

# Use max of both thresholds
# apply(X, 1, max) applies max() across each ROW (1 = rows, 2 = columns), so
# each system keeps whichever of its two thresholds is higher -- the binding
# constraint. na.rm=T ignores a missing tier rather than returning NA.
t <- apply(data.frame(cbind(t1, t2)), 1, max, na.rm = T)
# ignore warnings

# Replace -Inf with NA (when both values were NA)
# max() of nothing returns -Inf once the NAs have been removed, which is the
# warning referred to above. Those cases genuinely have no threshold recorded,
# so -Inf is converted back to NA.
t <- car::recode(t, "-Inf=NA")

realCases$highThreshold <- t # Ignore tiers 3 and 4 for the time being


## Recode seat allocation formulas to GBM-compatible labels
# Recode tier 1 formula
# car::recode() takes its instructions as a single text string, with each rule
# written 'old value'='new value' and rules separated by semicolons. Note how
# many recorded labels collapse onto one simulated category: every variant of
# the Hagenbach-Bischoff quota becomes 'hagenbachbischoff', and both 'Absolute
# Majority' and 'Alternative Vote' become 'abs majority'. Merging categories
# this way is the whole point -- the models know only the right-hand vocabulary.
# as.character() then ensures the result is plain text rather than a factor.
realCases$Formula <- as.character(car::recode(
  realCases$bg_tier1_formula,
  "'Absolute Majority'='abs majority';
'Alternative Vote'='abs majority';
'Block Vote'='plurality';
'dHondt'='dhondt';
'Droop quota'='droop';
'Droop quota with largest remainders'='droop';
'Hagenbach-Bischoff quota'='hagenbachbischoff';
'Hagenbach-Bischoff quota with highest average remainders'='hagenbachbischoff';
'Hagenbach-Bischoff quota with largest remainders'='hagenbachbischoff';
'Hare quota'='hare';
'Hare quota with largest remainders'='hare';
'Hare quota with highest average remainders'='hare';
'Imperiali quota'='imperiali';
'Imperialli quota with largest remainders'='imperiali';
'Limited Nomination'='plurality';
'Limited Vote'='plurality';
'Modified Sainte-Lague'='modsaintlague';
'Party Block Vote'='plurality';
'Reinforced Imperiali quota'='imperiali';
'Sainte-Lague'='saintelague';
'Single Nontransferable Vote'='plurality';
'Single Transferable Vote'='droop';
'Single-Member-District-Plurality'='plurality';
'Two-Round Absolute Majority'='abs majority';
'Two-Round Majority-Plurality'='abs majority';
'Two-Round Majority Runoff'  ='abs majority';
'Two-Round Qualified Majority'='abs majority';
'Two-Round Party Block Vote'='abs majority';
else=NA"
))

# Recode tier 2 formula
realCases$Formula2 <- as.character(car::recode(
  realCases$bg_tier2_formula,
  "'Hare quota with largest remainders'='hare';
'dHondt'='dhondt';
'Hare'='hare';
'Hare quota'='hare';
'Limited Nomination'='plurality';
'Limited Vote'='plurality';
'Droop quota'='droop';
'Droop quota with largest remainders'='droop';
'Hare quota with highest average remainders'='hare';
'Sainte-Lague'='saintelague';
'Hagenbach-Bischoff quota with largest remainders'='hagenbachbischoff';
'Two-Round Party Block Vote'='plurality';
'Hagenbach-Bischoff quota'='hagenbachbischoff';
'Reinforced Imperiali quota'='imperiali';
'Single Nontransferable Vote'='plurality';
'Single-Member-District-Plurality'='plurality';
'Two-Round Majority Runoff'='abs majority';
'Modified Sainte-Lague'='modsaintlague';
'Party Block Vote'='plurality';
else=NA"
))

## Other variables in realCases that we need to code
## None of these variables as such exist in realCases
# ballot_type: "closed", "flexible", "open"
# pool_level: "party_list", "party", "candidate"
# new.nvotes: "numeric"


## Recode ballot structure variables (ballot type, vote type, pooling level)
# Recode ballot type based on preferential voting label
realCases$ballot_type <- car::recode(
  realCases$preferential,
  " 'Open List'='open'
                                     ; 'Quasi List'='open'
                                     ; 'Latent List'='open'
                                     ; 'Flexible List'='flexible'
                                     ; 'Closed List'='closed'
                                     ; else=NA"
)

# Standardize ballot info label
realCases$ballotInfo <- car::recode(realCases$ballotInfo, " 'Dividual '='Dividual'
                                    ; 'DIvidual'='Dividual'")

# Recode number of votes field into simpler categories
realCases$votesN <- car::recode(
  realCases$votesN,
  "'Less than seats'='less seats';
                           'Total candidates'='total candidates';
                           'Total Candidates'='total candidates';
                           'Total seats'='total seats';
                           'Total Seats'='total seats';
                           'One '='one';
                           'One'='one';
                           'Two'='two';
                           else=NA"
)

# Reclassify fully "free" ballot structure: open ballot + individual votes + full number of votes
realCases$ballot_type <- ifelse(realCases$ballot_type == "open" & realCases$ballotInfo == "Dividual" & realCases$votesN == "total seats", "free", realCases$ballot_type)


## Recode pooling level and adjust tier-specific pooling info

# Assign general pooling level based on electoral rule type
realCases$pool_level <- car::recode(
  realCases$bg_elecrule,
  "'Two-Round System'='party';
                                    'List Proportional Representation'='party';
                                    'Mixed Independent'='party';
                                    'Alternative Vote'='party';
                                    'Single-Member-District-Plurality'='party';
                                    'Limited Nomination'='party';
                                    'Mixed Dependent'='party';
                                    'Block Vote'='candidate';
                                    'Party Block Vote'='party';
                                    'Limited Vote'='candidate';
                                    'Two-Round System'='party';
                                    'Mixed Independent '='party';
                                    'Single Transferable Vote'='candidate';
                                    'Single Nontransferable Vote'='candidate'; else=NA"
)


# New variables in "real cases" dataset: poolingTier1 and poolingTier2
# We use them to build alternative "pool_level" variables

# Recode pooling level by tier (alternative, more precise specification)
realCases$pool_level_tier1 <- car::recode(
  realCases$poolingTier1,
  "'Party'='party';
                                            'Candidate'='candidate';
                                            'Sub-Party'='subparty'; else=NA"
)


realCases$pool_level_tier2 <- car::recode(
  realCases$poolingTier2,
  "'Party'='party';
                                            'Candidate'='candidate';
                                            'Sub-Party'='subparty'; else=NA"
)

# There are 16 observations with pool_level_tier1=="subparty"
# We can incorporate these within pool_level

# Reclassify "subparty" pooling as "party_list" for tier 1
realCases$pool_level <- ifelse(realCases$pool_level_tier1 == "subparty", "party_list", realCases$pool_level)


########################################################################
# Average district magnitudes
########################################################################
# Non-compensatory mixed system: Size of chamber / Sum of # of districts = (tier1.seats + tier2.seats)/(tier1.dist + tier2.dist)
# For non-MM systems:  we use only the first tier for non-MM systems.
########################################################################


# Classify electoral systems into broad family categories
realCases$family <- as.character(car::recode(realCases$bg_elecrule, "'Two-Round System'='TRS';
                                    'List Proportional Representation'='LPR';
                                    'Alternative Vote'='AV';
                                    'Single-Member-District-Plurality'='SMDP';
                                    'Limited Nomination'='LN';
                                    'Mixed Dependent'='MxD';
                                    'Block Vote'='MNTV';
                                    'Party Block Vote'='PBV';
                                    'Limited Vote'='LV';
                                    'Two Round Majority Runoff'='TRS';
                                    'Mixed Independent'='MxI';
                                    'Mixed Independent '='MxI';
                                    'Single Transferable Vote'='STV';
                                    'Single Nontransferable Vote'='SNTV'; else='missing'"))

# Identify TR-MNTV cases: Two-Round systems with more than one vote per voter
realCases$family <- as.character(ifelse(realCases$family == "TRS" & realCases$votesN != "one", "TR-MNTV", realCases$family))

# Label "Fortified PR" as a separate family
realCases$family <- as.character(ifelse(realCases$bg_tier1_formula == "Fortified PR", "fort.PR", realCases$family))

# Create detailed family classification by combining formula family and ballot type
realCases$expanded.family <- realCases$family
realCases$expanded.family <- ifelse(realCases$ballot_type == "closed" &
  realCases$family == "LPR" &
  !is.na(realCases$family), "CLPR",
ifelse(realCases$ballot_type == "open" &
  realCases$family == "LPR" &
  !is.na(realCases$family), "OLPR",
ifelse(realCases$ballot_type == "flexible" &
  realCases$family == "LPR" &
  !is.na(realCases$family), "FLPR",
ifelse(realCases$ballot_type == "free" &
  realCases$family == "LPR" &
  !is.na(realCases$family), "FrLPR", realCases$family)
)
)
)

# Manually assign OLPR for countries where it is known but not well coded
realCases$expanded.family[is.na(realCases$expanded.family) & realCases$country %in% c("Colombia", "Iceland")] <- "OLPR"

# Drop cases without known family classification
## (to revise later)
realCases <- realCases[!is.na(realCases$expanded.family), ]

# Build variables to create Avg.M later
# Replace NA in upper chamber seats with 0
realCases$bg_upperseats <- car::recode(realCases$bg_upperseats, "NA=0")

# Derive total seats by tier
seats.1 <- ifelse(is.na(realCases$bg_seats), 0, realCases$bg_seats - realCases$bg_upperseats)
seats.2 <- ifelse(is.na(realCases$bg_upperseats), 0, realCases$bg_upperseats)
tier.1 <- ifelse(is.na(realCases$bg_tier1_districts), 0, realCases$bg_tier1_districts)
tier.2 <- ifelse(is.na(realCases$bg_tier2_districts), 0, realCases$bg_tier2_districts)
tier.3 <- ifelse(is.na(realCases$bg_tier3_districts), 0, realCases$bg_tier3_districts)
tier.4 <- ifelse(is.na(realCases$bg_tier4_districts), 0, realCases$bg_tier4_districts)

# # For mixed dependent systems, use the formula of the second tier (bg_tier2_formula).
realCases$Formula <- as.character(ifelse(realCases$family == "MxD" & !is.na(realCases$family),
  realCases$Formula2, realCases$Formula
))
realCases$Formula <- as.factor(realCases$Formula)

# Average M for non-mixed systems: seats / districts
avg.M <- seats.1 / tier.1

# Mixed independent: average over both tiers
avg.M.mixed.independent <- (seats.1 + seats.2) / (tier.1 + tier.2)
avg.M.mixed.independent <- ifelse(is.infinite(avg.M.mixed.independent), NA, avg.M.mixed.independent)
avg.M <- ifelse(realCases$family == "MxI" & !is.na(realCases$family), avg.M.mixed.independent, avg.M)

# Mixed dependent: only use upper tier (tier.2)
avg.M.mixed.dependent <- seats.2 / tier.2
avg.M.mixed.dependent <- ifelse(is.infinite(avg.M.mixed.dependent), NA, avg.M.mixed.dependent)
avg.M <- ifelse(realCases$family == "MxD" & !is.na(realCases$family), avg.M.mixed.dependent, avg.M)

# Fix infinite values and assign M = 1 for SMDP systems
avg.M <- ifelse(is.infinite(avg.M), NA, avg.M)
avg.M[realCases$family == "SMDP" & !is.na(realCases$family)] <- 1

# Store and overwrite M values
# Now M needs to get avg.M (we need to preserve "M" as a name)
# because the predictive functions below look for "M"
realCases$old.M <- realCases$M # Keep the originally-coded values of M elsewhere
realCases$M <- avg.M
realCases$avg.M <- avg.M


### Need a placeholder for M, for when we add a different set of calculations in the lines that follow

## Duplicate and split Mixed Independent systems into separate tiers

# Duplicate "mixed independent" (expanded.family=="MxI") observations
mixCases.tier1 <- realCases[realCases$expanded.family == "MxI", ]
mixCases.tier2 <- realCases[realCases$expanded.family == "MxI", ]
realCases <- realCases[-which(realCases$expanded.family == "MxI"), ]

# Note that this duplicates mixed independent observations
# This is resolved later

# Recalculate M for tier-specific observations
seats.tier1 <- ifelse(is.na(mixCases.tier1$bg_seats), 0, mixCases.tier1$bg_seats - mixCases.tier1$bg_upperseats)
tier.tier1 <- ifelse(is.na(mixCases.tier1$bg_tier1_districts), 0, mixCases.tier1$bg_tier1_districts)
seats.tier2 <- ifelse(is.na(mixCases.tier2$bg_upperseats), 0, mixCases.tier2$bg_upperseats)
tier.tier2 <- ifelse(is.na(mixCases.tier2$bg_tier2_districts), 0, mixCases.tier2$bg_tier2_districts)

mixCases.tier1$M <- ifelse(is.infinite(seats.tier1 / tier.tier1), NA, seats.tier1 / tier.tier1)
mixCases.tier2$M <- ifelse(is.infinite(seats.tier2 / tier.tier2), NA, seats.tier2 / tier.tier2)

mixCases.tier1$expanded.family <- "MxI.1"
mixCases.tier2$expanded.family <- "MxI.2"

# Check if three different datasets have identical column names (they should), then rbind them
identical(colnames(realCases), colnames(mixCases.tier1))

## [1] TRUE

identical(colnames(realCases), colnames(mixCases.tier2))

## [1] TRUE

# Combine split mixed-tier datasets with full dataset
realCases <- rbind(realCases, mixCases.tier1, mixCases.tier2)

# Final recoding for votes and numeric vote counts
realCases$votesN[realCases$expanded.family == "FLPR"] <- "one"
realCases$votesN[realCases$expanded.family == "OLPR"] <- "one"

# Tabulate vote types for reporting
table(realCases$votesN) # for descriptive table in Chapter 6

## 
##       less seats              one total candidates      total seats 
##                3             1509               48              131 
##              two 
##               22

# Translate votesN into numeric values for use in prediction
realCases$new.nvotes <- ifelse(!is.na(realCases$votesN) & realCases$votesN == "one", 1,
  ifelse(!is.na(realCases$votesN) & realCases$votesN == "two", 2,
    ifelse(!is.na(realCases$votesN) & realCases$votesN == "total candidates", realCases$M,
      ifelse(!is.na(realCases$votesN) & realCases$votesN == "total seats", realCases$M,
        ifelse(!is.na(realCases$votesN) & realCases$votesN == "less seats", realCases$M - 1, NA)
      )
    )
  )
)

# Fix known edge cases for Mixed + Dividual + Two vote systems (Correct realCases$new.nvotes)
mixed.elecrule <- grep("Mixed", realCases$bg_elecrule)
dividual <- grep("vidual", realCases$ballotInfo)
two.votesN <- grep("two", realCases$votesN)
mixedividual <- intersect(mixed.elecrule, dividual)
realCases$new.nvotes[intersect(mixedividual, two.votesN)] <- 1

# Set vote count to 1 in categorical ballots
realCases$new.nvotes[realCases$ballotInfo == "Categorical"] <- 1

## Final clean-up: drop invalid cases
# Drop Cases for which we don't have a valid value for M
# PCS: we still need to deal with MM conditional
realCases <- realCases[realCases$M != 0 & !is.na(realCases$M), ]
#  Drop cases with missing family classification (expanded.family is missing)
realCases <- realCases[realCases$expanded.family != "missing", ]

# Normalize threshold for GBM input (0–1 scale)
realCases$Threshold <- realCases$highThreshold / 100 # Threshold is in the unit-range in GBMs


###################################################
### Further corrections on real-cases database, ###
### to make comformable to GBM predictions      ###
###################################################


# Extract GBM models to retrieve factor level structures for compatibility checks
optimalGBMIntra <- totalAP.objects$optimalGBMIntra[[1]]
optimalGBM <- totalEffENP.objects$optimalGBM[[1]]

# PCS We don't have ballot_type == "free" in the sims
# Adjust "free" ballot types to "open" since GBM models do not include "free" as a valid level
realCases$ballot_type[realCases$ballot_type == "free"] <- "open"

### Further corrections to realCases to accommodate predictions
# Convert ballot_type to a factor with levels matching the GBM training model
realCases$ballot_type <- factor(realCases$ballot_type,
  levels = optimalGBM$var.levels[[4]]
)
# Convert pool_level to a factor with levels matching the GBM model
realCases$pool_level <- factor(realCases$pool_level,
  levels = optimalGBM$var.levels[[3]]
)

## Normalize system families to majoritarian, proportional, and mixed
# PCS: Create family to use in the predictions
realCases$family[realCases$family %in% c("MxD", "MxI")] <- "mixed"
realCases$family[realCases$family %in% c("fort.PR", "LPR", "STV")] <- "pr"
realCases$family[!(realCases$family %in% c("pr", "mixed"))] <- "maj"

# PCS: votes
# Create new.nvotes with cleaned categorical labels
realCases$new.nvotes <- realCases$votesN
realCases$new.nvotes[realCases$new.nvotes == "less seats"] <- "LessSeats"
realCases$new.nvotes[realCases$new.nvotes == "one"] <- "One"
realCases$new.nvotes[realCases$new.nvotes == "total candidates"] <- "TotalCandidates"
realCases$new.nvotes[realCases$new.nvotes == "total seats"] <- "TotalSeats"
# Special case: in TR-MNTV, treat 'two' as equivalent to 'TotalSeats'
realCases$new.nvotes[realCases$new.nvotes == "two" & realCases$expanded.family == "TR-MNTV"] <- "TotalSeats"
# In other systems, recode 'two' as 'One' (conservative assumption)
realCases$new.nvotes[realCases$new.nvotes == "two"] <- "One"
# Convert new.nvotes to factor with model-expected levels
realCases$new.nvotes <- factor(realCases$new.nvotes,
  levels = optimalGBM$var.levels[[2]]
)

##  Final predictor formatting for use in prediction
# Set formula and threshold variables using harmonized versions
realCases$formula <- realCases$Formula
realCases$threshold <- realCases$Threshold
# Replace NAs in threshold with 0 (default for countries without a threshold)
realCases$threshold <- car::recode(realCases$threshold, "NA=0")

The predictions of TDE and AP based on district- and country-level GBM predictive models come next.

## Using the country-level GBM object ##
totalEff.hat <- pers.hat <- list()

for (i in 1:5) {
  optimalGBMIntra <- totalAP.objects$optimalGBMIntra[[i]]
  optimalGBM <- totalEffENP.objects$optimalGBM[[i]]

  # Predict TDE and AP for all systems using country-level GBMs
  totalEff.hat[[i]] <- predict(optimalGBM, realCases,
    n.trees = optimalGBM$n.trees
  )

  pers.hat[[i]] <- predict(optimalGBMIntra, realCases,
    n.trees = optimalGBMIntra$n.trees
  )
}

# Store the five replicate predictions for TDE and AP
realCases$totalEff.hat1 <- totalEff.hat[[1]]
realCases$totalEff.hat2 <- totalEff.hat[[2]]
realCases$totalEff.hat3 <- totalEff.hat[[3]]
realCases$totalEff.hat4 <- totalEff.hat[[4]]
realCases$totalEff.hat5 <- totalEff.hat[[5]]

realCases$pers.hat1 <- pers.hat[[1]]
realCases$pers.hat2 <- pers.hat[[2]]
realCases$pers.hat3 <- pers.hat[[3]]
realCases$pers.hat4 <- pers.hat[[4]]
realCases$pers.hat5 <- pers.hat[[5]]

# Compute average across the five GBM replicates
realCases$totalEff.hat <- with(realCases, (totalEff.hat1 +
  totalEff.hat2 +
  totalEff.hat3 +
  totalEff.hat4 +
  totalEff.hat5) / 5)

realCases$pers.hat <- with(realCases, (pers.hat1 +
  pers.hat2 +
  pers.hat3 +
  pers.hat4 +
  pers.hat5) / 5)


## Using the district-level GBM object ##
totalEff.hat <- pers.hat <- list()

for (i in 1:5) {
  optimalGBMIntra <- totalAP.objects.district$optimalGBM[[i]]
  optimalGBM <- totalEffENP.objects.district$optimalGBM[[i]]

  totalEff.hat[[i]] <- predict(optimalGBM, realCases,
    n.trees = optimalGBM$n.trees
  )

  pers.hat[[i]] <- predict(optimalGBMIntra, realCases,
    n.trees = optimalGBMIntra$n.trees
  )
}

# Store district-level predictions
realCases$totalEff.hat1.district <- totalEff.hat[[1]]
realCases$totalEff.hat2.district <- totalEff.hat[[2]]
realCases$totalEff.hat3.district <- totalEff.hat[[3]]
realCases$totalEff.hat4.district <- totalEff.hat[[4]]
realCases$totalEff.hat5.district <- totalEff.hat[[5]]

realCases$pers.hat1.district <- pers.hat[[1]]
realCases$pers.hat2.district <- pers.hat[[2]]
realCases$pers.hat3.district <- pers.hat[[3]]
realCases$pers.hat4.district <- pers.hat[[4]]
realCases$pers.hat5.district <- pers.hat[[5]]

# Compute mean across the five district-based models
realCases$totalEff.hat.district <- with(realCases, (totalEff.hat1.district +
  totalEff.hat2.district +
  totalEff.hat3.district +
  totalEff.hat4.district +
  totalEff.hat5.district) / 5)

realCases$pers.hat.district <- with(realCases, (pers.hat1.district +
  pers.hat2.district +
  pers.hat3.district +
  pers.hat4.district +
  pers.hat5.district) / 5)


### For country-level GBM predictions, we need to combine M ###
## Adjust predictions for Mixed Independent (MxI) systems using weighted tiers
## Regroup MxI observations, eliminate double-observations

# Isolate MxI cases (duplicated in two rows per system) and create IDs
mixCases <- realCases[grep("MxI", realCases$expanded.family), ]
mixCases$id <- paste(mixCases$country, mixCases$year, sep = "-")


# Identify valid cases with both tiers present

# We need to drop 33 cases of MxI systems for which we lack district size for PR tiers
mixCases$id[mixCases$expanded.family == "MxI.1"][!is.element(mixCases$id[mixCases$expanded.family == "MxI.1"], mixCases$id[mixCases$expanded.family == "MxI.2"])]

##  [1] "Burkina Faso-1997" "Madagascar-1998"   "Madagascar-2002"  
##  [4] "Madagascar-2007"   "Madagascar-2013"   "Mauritania-1992"  
##  [7] "Mauritania-1996"   "Mauritania-2001"   "Morocco-1970"     
## [10] "Morocco-1977"      "Morocco-1984"      "Morocco-1993"     
## [13] "Nicaragua-1984"    "Nicaragua-1990"    "Panama-1984"      
## [16] "Panama-1989"       "Panama-1994"       "Panama-1999"      
## [19] "Panama-2004"       "Panama-2009"       "Panama-2014"

# But we can obtain average weights for 114 cases where we have the weights
changeOutcomeCountries <- mixCases$id[mixCases$expanded.family == "MxI.1"][is.element(mixCases$id[mixCases$expanded.family == "MxI.1"], mixCases$id[mixCases$expanded.family == "MxI.2"])]


# Combine predictions from both tiers using seat-based weights
totalEff.weighted <- c()
pers.weighted <- c()

for (i in 1:length(unique(mixCases$id))) {
  case <- unique(mixCases$id)[i]

  # Calculate weights based on seat shares
  seats1 <- mixCases$bg_seats[mixCases$id == case & mixCases$expanded.family == "MxI.1"] - mixCases$bg_upperseats[mixCases$id == case & mixCases$expanded.family == "MxI.1"] # unclear whehter whe should have bg_tier1_districts here, or bg_seats

  seats2 <- mixCases$bg_upperseats[mixCases$id == case & mixCases$expanded.family == "MxI.2"]
  weight1 <- seats1 / (seats1 + seats2)


  # ls.l <- mixCases$ls.lsq.hat[mixCases$id==case & mixCases$expanded.family=="MxI"]*weight1 +
  #   mixCases$ls.lsq.hat[mixCases$id==case & mixCases$expanded.family=="MxI.2"]*(1-weight1)
  # ls.p <- mixCases$ls.pers.hat[mixCases$id==case & mixCases$expanded.family=="MxI"]*weight1 +
  #   mixCases$ls.pers.hat[mixCases$id==case & mixCases$expanded.family=="MxI.2"]*(1-weight1)


  # Combine TDE and AP predictions using weights
  l <- mixCases$totalEff.hat[mixCases$id == case & mixCases$expanded.family == "MxI.1"] * weight1 +
    mixCases$totalEff.hat[mixCases$id == case & mixCases$expanded.family == "MxI.2"] * (1 - weight1)
  p <- mixCases$pers.hat[mixCases$id == case & mixCases$expanded.family == "MxI.1"] * weight1 +
    mixCases$pers.hat[mixCases$id == case & mixCases$expanded.family == "MxI.2"] * (1 - weight1)


  # l.sincere <- mixCases$lsq.hat.sincere[mixCases$id==case & mixCases$expanded.family=="MxI"]*weight1 +
  #              mixCases$lsq.hat.sincere[mixCases$id==case & mixCases$expanded.family=="MxI.2"]*(1-weight1)
  # p.sincere <- mixCases$pers.sincere[mixCases$id==case & mixCases$expanded.family=="MxI"]*weight1 +
  #              mixCases$pers.sincere[mixCases$id==case & mixCases$expanded.family=="MxI.2"]*(1-weight1)
  # l.strategic <- mixCases$lsq.hat.strategic[mixCases$id==case & mixCases$expanded.family=="MxI"]*weight1 +
  #                mixCases$lsq.hat.strategic[mixCases$id==case & mixCases$expanded.family=="MxI.2"]*(1-weight1)
  # p.strategic <- mixCases$pers.strategic[mixCases$id==case & mixCases$expanded.family=="MxI"]*weight1 +
  #                mixCases$pers.strategic[mixCases$id==case & mixCases$expanded.family=="MxI.2"]*(1-weight1)


  # Store weighted values
  totalEff.weighted <- c(l, totalEff.weighted)
  pers.weighted <- c(p, pers.weighted)


  # lsq.weighted  <- c(l, lsq.weighted)
  # pers.weighted <- c(p, pers.weighted)
  # lsq.sincere.weighted  <- c(l, lsq.sincere.weighted)
  # pers.sincere.weighted <- c(p, pers.sincere.weighted)
  # lsq.strategic.weighted  <- c(l, lsq.strategic.weighted)
  # pers.strategic.weighted <- c(p, pers.strategic.weighted)
}


## Replace MxI entries with weighted estimates

realCases$id <- paste(realCases$country, realCases$year, sep = "-")
changeOutcomeCountries[is.element(changeOutcomeCountries, realCases$id)]

##   [1] "Albania-1996"                      "Albania-1997"                     
##   [3] "Armenia-1995"                      "Armenia-1999"                     
##   [5] "Armenia-2003"                      "Armenia-2007"                     
##   [7] "Armenia-2012"                      "Azerbaijan-1995"                  
##   [9] "Azerbaijan-2000"                   "Bulgaria-1990"                    
##  [11] "Bulgaria-2009"                     "Cameroon-1992"                    
##  [13] "Cameroon-1997"                     "Cameroon-2002"                    
##  [15] "Cameroon-2007"                     "Cameroon-2013"                    
##  [17] "Chad-2002"                         "Chad-2011"                        
##  [19] "Croatia-1992"                      "Croatia-1995"                     
##  [21] "Democratic Republic of Congo-2006" "Democratic Republic of Congo-2011"
##  [23] "Egypt-1987"                        "Egypt-2011"                       
##  [25] "Georgia-1992"                      "Georgia-1995"                     
##  [27] "Georgia-1999"                      "Georgia-2003"                     
##  [29] "Georgia-2004"                      "Georgia-2008"                     
##  [31] "Georgia-2012"                      "Guinea-1995"                      
##  [33] "Guinea-2002"                       "Guinea-2013"                      
##  [35] "Iceland-1946"                      "Iceland-1949"                     
##  [37] "Iceland-1953"                      "Iceland-1956"                     
##  [39] "Iceland-1959"                      "Japan-1996"                       
##  [41] "Japan-2000"                        "Japan-2003"                       
##  [43] "Japan-2005"                        "Japan-2009"                       
##  [45] "Japan-2012"                        "Japan-2014"                       
##  [47] "Jordan-2013"                       "Kazakhstan-1999"                  
##  [49] "Kazakhstan-2004"                   "Kyrgyzstan-2000"                  
##  [51] "Lithuania-1992"                    "Lithuania-1996"                   
##  [53] "Lithuania-2000"                    "Lithuania-2004"                   
##  [55] "Lithuania-2008"                    "Lithuania-2012"                   
##  [57] "Macedonia-1998"                    "Mexico-1979"                      
##  [59] "Mexico-1982"                       "Mexico-1985"                      
##  [61] "Mexico-1988"                       "Mexico-1991"                      
##  [63] "Mexico-1994"                       "Mexico-1997"                      
##  [65] "Mexico-2000"                       "Mexico-2003"                      
##  [67] "Mexico-2006"                       "Mexico-2009"                      
##  [69] "Mexico-2012"                       "Mexico-2015"                      
##  [71] "Mongolia-2012"                     "Nepal-2008"                       
##  [73] "Nepal-2013"                        "Nicaragua-1996"                   
##  [75] "Nicaragua-2001"                    "Nicaragua-2006"                   
##  [77] "Nicaragua-2011"                    "Niger-1993"                       
##  [79] "Niger-1995"                        "Niger-1996"                       
##  [81] "Niger-1999"                        "Niger-2004"                       
##  [83] "Niger-2009"                        "Niger-2011"                       
##  [85] "Russia-1993"                       "Russia-1995"                      
##  [87] "Russia-1999"                       "Russia-2003"                      
##  [89] "Senegal-1983"                      "Senegal-1988"                     
##  [91] "Senegal-1993"                      "Senegal-1998"                     
##  [93] "Senegal-2001"                      "Senegal-2007"                     
##  [95] "Senegal-2012"                      "Serbia-1992"                      
##  [97] "South Korea-1963"                  "South Korea-1967"                 
##  [99] "South Korea-1971"                  "South Korea-1981"                 
## [101] "South Korea-1985"                  "South Korea-1988"                 
## [103] "South Korea-1992"                  "South Korea-1996"                 
## [105] "South Korea-2000"                  "South Korea-2004"                 
## [107] "South Korea-2008"                  "South Korea-2012"                 
## [109] "Sudan-2010"                        "Sudan-2015"                       
## [111] "Taiwan-1992"                       "Taiwan-1995"                      
## [113] "Taiwan-1998"                       "Taiwan-2001"                      
## [115] "Taiwan-2004"                       "Taiwan-2008"                      
## [117] "Taiwan-2012"                       "Tajikistan-2000"                  
## [119] "Tajikistan-2005"                   "Tajikistan-2010"                  
## [121] "Tajikistan-2015"                   "Thailand-2001"                    
## [123] "Thailand-2005"                     "Thailand-2006"                    
## [125] "Thailand-2007"                     "Thailand-2011"                    
## [127] "Thailand-2014"                     "Ukraine-1998"                     
## [129] "Ukraine-2002"                      "Ukraine-2012"                     
## [131] "Ukraine-2014"                      "Venezuela-2010"                   
## [133] "Venezuela-2015"                    "Zimbabwe-2013"

# Use this function to find elements in realCases that are MxI,
# therefore have a weighted totalEff and weighted pers score in
# totalEff.weighted and pers.weighted

# Match weighted results to MxI entries in realCases
findCorrespondence <- function(name.a, list.b) {
  pos <- which(list.b == name.a)
  return(pos)
}

# Take the appropriate totalEff.weighted and pers.weighted scores and
# substitute them in MxI systems in realCases

# Apply the weighted scores
for (i in 1:nrow(realCases)) {
  case <- realCases$id[i]
  if (!is.element(case, changeOutcomeCountries)) {
    next
  } else {
    realCases$totalEff.hat[i] <- totalEff.weighted[findCorrespondence(name.a = case, list.b = changeOutcomeCountries)]
    realCases$pers.hat[i] <- pers.weighted[findCorrespondence(name.a = case, list.b = changeOutcomeCountries)]
  }
}

# Remove "spare" MxI.2 entries; restore expanded.family==MxI.1 to expanded.family==MxI

# Remove duplicated tier-2 entries (MxI.2) and standardize naming
realCases <- realCases[!is.element(realCases$expanded.family, "MxI.2"), ]

# Recode expanded family for display/plotting purposes
realCases$expanded.family <- car::recode(realCases$expanded.family, "'MxI.1'='MMI'
                                          ; 'MxD'='MMC'
                                          ; 'TRS'='TR'")

The following chunk just adds some useful notation (for example, country name abbreviations).

## Add country codes and abbreviations for plotting

## Add 2-letter ISO codes to each observation
realCases$iso2 <- car::recode(realCases$country, "'Afghanistan'='AF';
                              'Angola'='AO';
                              'Austria'='AT';
                              'Belarus'='BY';
                              'Bhutan'='BT';
                              'Brazil'='BR';
                              'Burma/Myanmar'='MM';
                              'Canada'='CA';
                              'Chile'='CL';
                              'Croatia'='HR';
                              'Czech Republic'='CZ';
                              'Dominican Republic'='DO';
                              'El Salvador'='SV';
                              'Finland'='FI';
                              'German Democratic Republic'='DD';
                              'Greece'='GR';
                              'Haiti'='HT';
                              'Iceland'='IS';
                              'Iran'='IR';
                              'Israel'='IL';
                              'Japan'='JP';
                              'Kenya'='KE';
                              'Lebanon'='LB';
                              'Libya'='LY';
                              'Madagascar'='MG';
                              'Mali'='ML';
                              'Moldova'='MD';
                              'Morocco'='MA';
                              'Nepal'='NP';
                              'Nicaragua'='NI';
                              'Oman'='OM';
                              'Paraguay'='PY';
                              'Poland'='PL';
                              'Romania'='RO';
                              'Senegal'='SN';
                              'Singapore'='SG';
                              'Somalia'='SO';
                              'Spain'='ES';
                              'Swaziland'='SZ';
                              'Syria'='SY';
                              'Tanzania'='TZ';
                              'Togo'='TG';
                              'Turkmenistan'='TM';
                              'United Kingdom'='UK';
                              'Venezuela'='VE';
                              'Zimbabwe'='ZW';
                              'Georgia'='GE';
                              'Mauritania'='MR';
                              'Albania'='AL';
                              'Argentina'='AR';
                              'Azerbaijan'='AZ';
                              'Belgium'='BE';
                              'Bolivia'='BO';
                              'Bulgaria'='BG';
                              'Burundi'='BI';
                              'Central African Republic'='CF';
                              'Colombia'='CO';
                              'Cuba'='CU';
                              'Democratic Republic of Congo'='CD';
                              'Ecuador'='EC';
                              'Estonia'='EE';
                              'France'='FR';
                              'Germany'='DE';
                              'Guatemala'='GT';
                              'Honduras'='HN';
                              'India'='IN';
                              'Iraq'='IQ';
                              'Italy'='IT';
                              'Jordan'='JO';
                              'Kyrgyzstan'='KG';
                              'Lesotho'='LS';
                              'Luxembourg'='LU';
                              'Malawi'='MW';
                              'Mauritius'='MU';
                              'Mongolia'='MN';
                              'Mozambique'='MZ';
                              'Netherlands'='NL';
                              'Nigeria'='NG';
                              'Pakistan'='PK';
                              'Peru'='PE';
                              'Portugal'='PT';
                              'Russia'='RU';
                              'Serbia'='RS';
                              'Slovakia'='SK';
                              'South Africa'='ZA';
                              'Sri Lanka'='LK';
                              'Sweden'='SE';
                              'Taiwan'='TW';
                              'Thailand'='TH';
                              'Tunisia'='TN';
                              'Uganda'='UG';
                              'United States of America'='US';
                              'Yemen'='YE';
                              'Armenia'='AM';
                              'Guinea'='GN';
                              'Niger'='NE';
                              'Algeria'='DZ';
                              'Australia'='AU';
                              'Bangladesh'='BD';
                              'Benin'='BJ';
                              'Botswana'='BW';
                              'Burkina Faso'='BF';
                              'Cambodia'='KH';
                              'Chad'='TD';
                              'Costa Rica'='CR';
                              'Cyprus'='CY';
                              'Denmark'='DK';
                              'Egypt'='EG';
                              'Ethiopia'='ET';
                              'Gabon'='GA';
                              'Ghana'='GH';
                              'Guinea-Bissau'='GW';
                              'Hungary'='HU';
                              'Indonesia'='ID';
                              'Ireland'='IE';
                              'Ivory Coast'='CI';
                              'Kazakhstan'='KZ';
                              'Latvia'='LV';
                              'Liberia'='LR';
                              'Macedonia'='MK';
                              'Malaysia'='MY';
                              'Mexico'='MX';
                              'Montenegro'='ME';
                              'Namibia'='NA';
                              'New Zealand'='NZ';
                              'Norway'='NO';
                              'Panama'='PA';
                              'Philippines'='PH';
                              'Republic of the Congo'='CG';
                              'Rwanda'='RW';
                              'Sierra Leone'='SL';
                              'Slovenia'='SI';
                              'South Korea'='KR';
                              'Sudan'='SD';
                              'Switzerland'='CH';
                              'Tajikistan'='TJ';
                              'The Gambia'='GM';
                              'Turkey'='TR';
                              'Ukraine'='UA';
                              'Uruguay'='UY';
                              'Zambia'='ZM';
                              'Cameroon'='CM';
                              'Republic of Vietnam'='VN';
                              'Lithuania'='LT'")


# Add 3-letter country codes for alternative labeling
realCases$scode <- car::recode(realCases$country, "'Afghanistan'='AFG';
                              'Angola'='ANG';
                              'Austria'='AUS';
                              'Belarus'='BLR';
                              'Bhutan'='BHU';
                              'Brazil'='BRA';
                              'Burma/Myanmar'='MYA';
                              'Canada'='CAN';
                              'Chile'='CHL';
                              'Croatia'='CRO';
                              'Czech Republic'='CZR';
                              'Dominican Republic'='DOM';
                              'El Salvador'='SAL';
                              'Finland'='FIN';
                              'German Democratic Republic'='GDR';
                              'Greece'='GRC';
                              'Haiti'='HAI';
                              'Iceland'='ICE';
                              'Iran'='IRN';
                              'Israel'='ISR';
                              'Japan'='JPN';
                              'Kenya'='KEN';
                              'Lebanon'='LEB';
                              'Libya'='LIB';
                              'Madagascar'='MAG';
                              'Mali'='MLI';
                              'Moldova'='MLD';
                              'Morocco'='MOR';
                              'Nepal'='NEP';
                              'Nicaragua'='NIC';
                              'Oman'='OMA';
                              'Paraguay'='PAR';
                              'Poland'='POL';
                              'Romania'='RUM';
                              'Senegal'='SEN';
                              'Singapore'='SIN';
                              'Somalia'='SOM';
                              'Spain'='SPN';
                              'Swaziland'='SWA';
                              'Syria'='SYR';
                              'Tanzania'='TAZ';
                              'Togo'='TOG';
                              'Turkmenistan'='TKM';
                              'United Kingdom'='UKG';
                              'Venezuela'='VEN';
                              'Zimbabwe'='ZIM';
                              'Georgia'='GRG';
                              'Mauritania'='MAA';
                              'Albania'='ALB';
                              'Argentina'='ARG';
                              'Azerbaijan'='AZE';
                              'Belgium'='BEL';
                              'Bolivia'='BOL';
                              'Bulgaria'='BUL';
                              'Burundi'='BUI';
                              'Central African Republic'='CEN';
                              'Colombia'='COL';
                              'Cuba'='CUB';
                              'Democratic Republic of Congo'='ZAI';
                              'Ecuador'='ECU';
                              'Estonia'='EST';
                              'France'='FRN';
                              'Germany'='GMY';
                              'Guatemala'='GUA';
                              'Honduras'='HON';
                              'India'='IND';
                              'Iraq'='IRQ';
                              'Italy'='ITA';
                              'Jordan'='JOR';
                              'Kyrgyzstan'='KYR';
                              'Lesotho'='LES';
                              'Luxembourg'='LUX';
                              'Malawi'='MAW';
                              'Mauritius'='MAS';
                              'Mongolia'='MON';
                              'Mozambique'='MZM';
                              'Netherlands'='NTH';
                              'Nigeria'='NIG';
                              'Pakistan'='PAK';
                              'Peru'='PER';
                              'Portugal'='POR';
                              'Russia'='RUS';
                              'Serbia'='SER';
                              'Slovakia'='SLO';
                              'South Africa'='SAF';
                              'Sri Lanka'='SRI';
                              'Sweden'='SWD';
                              'Taiwan'='TAW';
                              'Thailand'='THI';
                              'Tunisia'='TUN';
                              'Uganda'='UGA';
                              'United States of America'='USA';
                              'Yemen'='YEM';
                              'Armenia'='ARM';
                              'Guinea'='GUI';
                              'Niger'='NIR';
                              'Algeria'='ALG';
                              'Australia'='AUL';
                              'Bangladesh'='BNG';
                              'Benin'='BEN';
                              'Botswana'='BOT';
                              'Burkina Faso'='BFO';
                              'Cambodia'='CAM';
                              'Chad'='CHA';
                              'Costa Rica'='COS';
                              'Cyprus'='CYP';
                              'Denmark'='DEN';
                              'Egypt'='EGY';
                              'Ethiopia'='ETH';
                              'Gabon'='GAB';
                              'Ghana'='GHS';
                              'Guinea-Bissau'='GNB';
                              'Hungary'='HUN';
                              'Indonesia'='INS';
                              'Ireland'='IRE';
                              'Ivory Coast'='IVO';
                              'Kazakhstan'='KZK';
                              'Latvia'='LAT';
                              'Liberia'='LBR';
                              'Macedonia'='MAC';
                              'Malaysia'='MAL';
                              'Mexico'='MEX';
                              'Montenegro'='MNT';
                              'Namibia'='NAM';
                              'New Zealand'='NEW';
                              'Norway'='NOR';
                              'Panama'='PAN';
                              'Philippines'='PHI';
                              'Republic of the Congo'='CON';
                              'Rwanda'='RWA';
                              'Sierra Leone'='SIE';
                              'Slovenia'='SLV';
                              'South Korea'='KOR';
                              'Sudan'='SUD';
                              'Switzerland'='SWZ';
                              'Tajikistan'='TAJ';
                              'The Gambia'='GAM';
                              'Turkey'='TUR';
                              'Ukraine'='UKR';
                              'Uruguay'='URU';
                              'Zambia'='ZAM';
                              'Cameroon'='CAO';
                              'Republic of Vietnam'='VIE';
                              'Lithuania'='LIT'")


# Create a compact ID using iso2 code and year suffix
realCases$id2 <- paste0(realCases$iso2, substr(realCases$year, 3, 4))

# Create full ID with 3-letter code and full year (e.g., FRA1986)
realCases$full.id2 <- paste0(realCases$scode, realCases$year)


# Create a list of descriptive names for electoral family categories used in plots
familyTitle <- c(
  "Single Non-Transferable Vote",
  "Two-Round Runoff",
  "Mixed-Member Compensatory",
  "Closed-List Proportional Representation",
  "Limited Nomination",
  "Single-Member District Plurality",
  "Alternative Vote",
  "Flexible-List Proportional Representation",
  "Fortified PR",
  "Open-List Proportional Representation",
  "Multiple Non-Transferable Vote",
  "Free-List Proportional Representation",
  "Two-Round Multiple Non-Transferable Vote",
  "Party Block Vote",
  "Single Transferable Vote",
  "Limited Vote",
  "Mixed-Member Independent"
)

The following chunk adds information about the democracy/autocracy status of different country/year observations.

Requires PolityIV.RData. This block reads a cached extract of the Polity IV dataset that must sit in the same folder as this .Rmd file. The script make-PolityIV-cache.R, distributed alongside these files, downloads the data and writes that cache; it needs to be run only once.

The download was previously performed here, inline, by democracyData::download_polity_annual(). It was moved out for two reasons. First, it was the only step in all ten chapters that reached out to the internet at knit time, so the chapter could not be reproduced offline. Second, and more seriously, that function no longer returns the data the book was built on: the Polity IV annual series has been archived upstream and the function now downloads Polity 5, which revised its scores. Re-running the old code today could therefore reclassify country-years and silently change which cases appear in Figure 6.2.

# Install and load the 'democracyData' package for Polity IV data
# remotes::install_github("xmarquez/democracyData")

# The democracyData package is no longer loaded here. It is needed only to BUILD
# the cache, which is done once by make-PolityIV-cache.R; knitting this chapter
# needs nothing but the saved file.
# library (democracyData)

# Load the cached Polity IV annual dataset
# load() restores the object saved by make-PolityIV-cache.R under its original
# name, `Polity`, so every line below is unchanged from the original code.
# No path is given because knitr runs each chunk with the working directory set
# to the folder holding this .Rmd file, which is the repository root.
load("data/shared/PolityIV.RData")
# Create a unique country-year identifier to match with realCases
# Pasting country code and year together makes a single key that uniquely
# identifies an observation in both datasets, which is what the join needs.
Polity$full.id2 <- paste0(Polity$scode, Polity$year)
# Keep only relevant columns (ID and polity2 score)
Polity <- Polity %>% dplyr::select(full.id2, polity2)

# Merge democracy data into the main dataset
# left_join() keeps every row of realCases and attaches the matching Polity
# score where one exists, leaving NA where none does. "Left" is the important
# word: no election is dropped for lacking a Polity score.
realCases <- left_join(realCases, Polity)
# Define democracy dummy: Polity2 score > 6
# ifelse(test, yes, no) applied to the whole column at once: 1 for democracies,
# 0 otherwise. The threshold of 6 is the conventional Polity cut-point.
realCases$democ <- ifelse(realCases$polity2 > 6, 1, 0)
# Force democratic classification for well-established cases with missing Polity values
# These four are uncontroversially democratic over the period but have gaps in
# the Polity series; without this line they would be silently excluded from
# Figure 6.2 by the NA rather than by any substantive judgement.
realCases$democ[realCases$country %in% c("Germany", "Iceland", "Japan", "South Korea")] <- 1

The next chunk provides some ancillary information that will appear in plots.

From 1,528 elections to seventeen points. The figures do not plot every election. They plot one point per electoral system family — CLPR, OLPR, STV, SNTV and so on — placed at the average TDE and AP of all the real elections that used it, and surrounded by an ellipse showing how much those elections varied. This block computes those averages, spreads and counts.

# Aggregate predicted TDE and AP scores by electoral family (based on country-level GBM models)

#### Build I-I plot with real cases ####
# Build dataset averaging at the "electoral system" level,
# using the country-level GBM predictor object
# group_by() + summarize() collapses the data to one row per family, carrying
# four numbers: the mean and standard deviation on each of the two dimensions.
# The means fix where each family's point sits; the standard deviations set the
# size of the ellipse drawn around it. na.rm=T ignores missing scores rather
# than returning NA for the whole family.
mean_realCases <- realCases %>%
  group_by(expanded.family) %>%
  summarize(
    mean.totalEff.hat = mean(totalEff.hat, na.rm = T),
    sd.totalEff.hat = sd(totalEff.hat, na.rm = T),
    mean.pers.hat = mean(pers.hat, na.rm = T),
    sd.pers.hat = sd(pers.hat, na.rm = T)
  )

# Drop cases with undefined family classification
# Systems that could not be classified form a "missing" group, which is not a
# real family and so is excluded from the figures.
mean_realCases <- mean_realCases[mean_realCases$expanded.family != "missing", ]
# Count the number of system instances in each electoral family
# table() counts elections per family. These counts are reported in the figures
# so a reader can see which points rest on many cases and which on few.
instancesFamily <- table(realCases$expanded.family)
# Add instance counts to the aggregated dataset
# The bracketed lookup reorders the counts to match the row order of
# mean_realCases before cbind() attaches them as a new column.
mean_realCases <- cbind(instancesFamily[mean_realCases$expanded.family], mean_realCases)


# Define color scheme for main I-I plots (used in Figures 6.1–6.4)
# One colour per electoral family, in the same order as the rows of
# mean_realCases. Only greys are used, because the book prints in black and
# white; the commented-out scheme further below is a colour version kept for
# on-screen work.
color.scheme <- c(
  "black", "darkslategray", "darkslategray", "slategray", "black",
  "darkslategray", "black", "slategray", "slategray",
  "black", "darkslategray", "black", "black", "black",
  "darkslategray", "black", "black"
)

# Alternate color scheme for another plotting variant
color.scheme.2 <- c(
  "black", "black", "black", "black", "black",
  "black", "slategray", "slategray", "black",
  "black", "black", "black", "black", "black",
  "black", "black", "black"
)


# color.scheme <- c("black","red","red","purple","black"
#                   ,"blue","blue","purple","orange","orange"
#                   ,"purple","blue","purple","blue","red")

# Create lowercase acronyms for electoral families for use in annotations or legends
expanded.family.acronyms <- mean_realCases$expanded.family
expanded.family.acronyms <- car::recode(
  expanded.family.acronyms,
  "'CLPR'='clpr';
                                          'FLPR'='flpr';
                                          'fort.PR'='fort.pr';
                                          'FrLPR'='frlpr';
                                          'OLPR'='olpr';
                                          'STV'='stv'"
)

Finally, we draw the plots that appear throughout Chapter 6.

Figure 6.1: Model-Based Predicted TDE and AP Values

Figure 6.1 maps average electoral system families within the Interparty–Intraparty (I–I) space using predictions from gradient boosting machine (GBM) models. These plots visualize how different electoral systems shape incentives along two key dimensions: the Total Duvergerian Effect (TDE), representing interparty competition, and personalism scores, capturing intraparty competition.

The following chunk produces Figure 6.1b, which plots average predicted values for each electoral system family in the Interparty–Intraparty (I–I) space using country-level GBM models. It visualizes how electoral rules shape interparty and intraparty incentives at the national level.

Base R graphics, built up in layers. These figures use R’s built-in plotting system rather than ggplot2. It works like painting on a canvas: plot() opens the canvas, and each later command adds something to it. Nothing appears unless it is explicitly drawn, and the order of commands is the order of the layers.

The sequence here is worth following, because all five figures in this chapter repeat it:

plot(..., type="n", axes=F) — open a canvas of the right size but draw nothing. type="n" means “no points”; axes=F suppresses the default axes. This reserves the space so everything else can be placed by hand.
mtext() — write the axis labels in the margins.
axis(1) / axis(2) — draw the bottom and left axes (1 = bottom, 2 = left).
text() — write the grey quadrant descriptors inside the plot area.
points(ellipse(...)) — draw one dispersion ellipse per family.
text() again — finally, the family acronyms on top.

par(mar = c(4,4,1,2)) sets the margins in lines of text, always in the order bottom, left, top, right.

par(mar = c(4, 4, 1, 2))

# Initializes the base plot with no points (type="n") and custom axes/margins.
# Uses the mean values of predicted TDE and AP across electoral system families based on country-level GBM models.

# with(data, plot(...)) lets the columns be named directly. xlim and ylim fix
# the extent of the two dimensions so that Figures 6.1a and 6.1b share a common
# frame and can be compared against each other.
# with (mean_realCases, plot (x=mean.lsq.hat, y=mean.pers.hat
with(mean_realCases, plot(
  x = mean.totalEff.hat, # X-axis is the interparty dimension (TDE)
  y = mean.pers.hat # Y-axis is the intraparty dimension (AP)
  , xlab = "",
  ylab = "",
  xlim = c(0, 4), ylim = c(0.3, 0.6),
  axes = F, type = "n"
))


# Adds axis labels to convey the conceptual meaning of both dimensions.
mtext("Interparty dimension", side = 1, line = 2)
mtext("(Total Duvergerian effect)", side = 1, line = 3, cex = 0.8)
# mtext ("(average weighted variance)", side=1, line=3, cex=0.8)
# mtext ("(average disproportionality score)", side=1, line=3, cex=0.8)
mtext("Intraparty dimension", side = 2, line = 3)
mtext("(Personalism score)", side = 2, line = 2, cex = 0.8)
# Draws the default x and y axes.
axis(1)
axis(2)
# Adds quadrant descriptors for the I–I space:
# Y-axis describes levels of personalism
# X-axis describes how constraining the electoral system is on interparty competition
text(
  x = rep(2, 4), y = c(0.33, 0.41, 0.49, 0.57),
  labels = c("Highly Centralized", "Party-Oriented", "Candidate-Oriented", "Highly Individualistic"),
  col = "gray"
)
text(
  x = c(0.5, 1.75, 2.75, 3.5), y = rep(0.45, 4),
  labels = c("Very Weak", "Permissive", "Constraining", "Very Strong"),
  col = "gray", srt = 90
)
# Draws a dashed ellipse around each family’s mean point, with radius equal to half a standard deviation in both dimensions.
# AV and SNTV are excluded due to insufficient dispersion or atypical values.
# The ellipse is a visual summary of spread, not a confidence region: it is
# centred on the family mean and extends half a standard deviation on each
# dimension. A wide ellipse means elections using that family landed in very
# different places in the I–I space; a tight one means the family reliably
# produces the same incentives.
# is.element(a, b) asks whether a appears in b -- the same test as %in%.
# `next` skips to the next pass of the loop without drawing anything, which is
# how AV and SNTV are left out.
for (i in 1:nrow(mean_realCases)) {
  if (is.element(mean_realCases$expanded.family[i], c("AV", "SNTV"))) {
    next
  } else {
    points(
      ellipse(
        x = 0,
        centre = c(mean_realCases$mean.totalEff.hat[i], mean_realCases$mean.pers.hat[i]),
        scale = c(mean_realCases$sd.totalEff.hat[i] / 2, mean_realCases$sd.pers.hat[i] / 2)
      ),
      col = "gray", lwd = 1, lty = 2, type = "l"
    )
  }
}
# Labels each electoral family at its mean location using short acronyms and a predefined color scheme (majoritarian, mixed, proportional).
with(mean_realCases, text(xy.coords(mean.totalEff.hat, mean.pers.hat), labels = expanded.family.acronyms, col = color.scheme.2))
# Adds a legend to identify the electoral system families by color.
legend("bottomright",
  pch = 19, col = c("black", "slategray", "black"),
  legend = c("MAJORITARIAN", "Mixed", "Proportional"), bty = "n"
)

This chunk generates Figure 6.1a, based on predictions from district-level GBM models. It displays how the same electoral system families are positioned in the I–I space when accounting for subnational variation in institutional design.

# Aggregates district-level predictions (TDE and AP) by electoral family to produce average coordinates for each system type.

mean_realCases_district <- realCases %>%
  group_by(expanded.family) %>%
  summarize(
    mean.totalEff.hat = mean(totalEff.hat.district, na.rm = T),
    sd.totalEff.hat = sd(totalEff.hat.district, na.rm = T),
    mean.pers.hat = mean(pers.hat.district, na.rm = T),
    sd.pers.hat = sd(pers.hat.district, na.rm = T)
  )
mean_realCases_district <- mean_realCases_district[mean_realCases_district$expanded.family != "missing", ]

# Removes missing family classifications and appends the instance count from earlier preprocessing.
mean_realCases_district <- cbind(instancesFamily[mean_realCases_district$expanded.family], mean_realCases_district)

# Creates the same plotting space as Figure 6.1b, but this time using district-level GBM projections.
par(mar = c(4, 4, 1, 2))
# with (mean_realCases, plot (x=mean.lsq.hat, y=mean.pers.hat
with(mean_realCases_district, plot(
  x = mean.totalEff.hat, y = mean.pers.hat,
  xlab = "",
  ylab = "",
  xlim = c(0, 4), ylim = c(0.3, 0.6),
  axes = F, type = "n"
))
# Repeats axes and labels for consistent presentation.
mtext("Interparty dimension", side = 1, line = 2)
mtext("(Total Duvergerian effect)", side = 1, line = 3, cex = 0.8)
# mtext ("(average weighted variance)", side=1, line=3, cex=0.8)
# mtext ("(average disproportionality score)", side=1, line=3, cex=0.8)
mtext("Intraparty dimension", side = 2, line = 3)
mtext("(Personalism score)", side = 2, line = 2, cex = 0.8)
axis(1)
axis(2)
# Adds quadrant labels as in 6.1b, slightly shifting the x-coordinates for visual balance.
text(
  x = rep(2, 4), y = c(0.33, 0.41, 0.49, 0.57),
  labels = c("Highly Centralized", "Party-Oriented", "Candidate-Oriented", "Highly Individualistic"),
  col = "gray"
)
text(
  x = c(0.5, 1.25, 2.75, 3.5), y = rep(0.45, 4),
  labels = c("Very Weak", "Permissive", "Constraining", "Very Strong"),
  col = "gray", srt = 90
)
# Draws ellipses representing variation around each family’s average I–I coordinates, again skipping AV and SNTV.
for (i in 1:nrow(mean_realCases_district)) {
  if (is.element(mean_realCases_district$expanded.family[i], c("AV", "SNTV"))) {
    next
  } else {
    points(
      ellipse(
        x = 0,
        centre = c(mean_realCases_district$mean.totalEff.hat[i], mean_realCases_district$mean.pers.hat[i]),
        scale = c(mean_realCases_district$sd.totalEff.hat[i] / 2, mean_realCases_district$sd.pers.hat[i] / 2)
      ),
      col = "gray", lwd = 1, lty = 2, type = "l"
    )
  }
}
# Labels each electoral system type at its district-level mean position.
with(mean_realCases_district, text(xy.coords(mean.totalEff.hat, mean.pers.hat), labels = expanded.family.acronyms, col = color.scheme.2))
# Adds the same legend for comparison with the country-level plot.
legend("bottomright",
  pch = 19, col = c("black", "slategray", "black"),
  legend = c("MAJORITARIAN", "Mixed", "Proportional"), bty = "n"
)

These plots display the location of real-world electoral systems in the interparty–intraparty (I–I) space, as predicted by our simulation-based GBM models. The x-axis shows the Total Duvergerian Effect (TDE), our measure of interparty incentives, while the y-axis captures the personalism score, which reflects intraparty incentives. Figure 6.1b uses the country-level GBM predictions, while Figure 6.1a is based on the district-level model. Each electoral system family is represented by its average predicted values, with dashed ellipses indicating the standard deviation around those means. As discussed in the chapter, systems vary substantially in how they structure both dimensions of political competition, with proportional systems tending toward low interparty constraint and majoritarian systems clustering on the right side of the x-axis.

These plots highlight the multidimensional nature of electoral incentives and demonstrate how institutional rules combine to produce diverse political environments across the world.

How to read Figure 6.1. The plot is a map, not a graph of one variable against another: horizontal position is how strongly the system constrains competition between parties, vertical position how strongly it rewards individual candidates within them. The grey words are signposts for the regions, not data.

Each acronym marks one electoral family, sitting at its average across every real election that used it, with a dashed ellipse showing the spread. Three things reward attention: which quadrant a family occupies; how large its ellipse is, since a wide one means the same family produced very different incentives in different countries; and whether families overlap, which would mean the family label alone does not pin down the incentives.

Panels (a) and (b) show the same families computed two ways — district-level models in (a), country-level in (b). Families that move between the two panels are those where the aggregation choice of Table 6.2 matters most.

Ordering. Figure 6.1 is generated as panel (b) first and panel (a) second, so the code blocks appear in reverse order relative to the book. Figure 6.2 is produced last in this file, after Figures 6.3, 6.4 and 6.5.

Figure 6.3: Select Majoritarian Electoral Systems

Figure 6.3 compares a hand-picked set of majoritarian electoral systems within the Interparty–Intraparty (I–I) space, illustrating how electoral incentives vary across time and context even within this single family of systems. Each system appears twice: once based on the country-level GBM model (bold italic labels) and once based on the district-level GBM model (regular font). This dual representation shows that despite certain systems—such as single-member district plurality (SMDP)—receiving identical simulation scores across models, the broader majoritarian family generates substantial variation in interparty and intraparty incentives. Note that this plot (and also plots in Figures 6.4 and 6.5) may be slightly different from the one in the book because we employ a random “jitter” statement on our code that adds some noise to point coordinates so that the election labels can be read more clearly. Beyond label readability, this change is inconsequential.

Individual elections, not families. Figures 6.3 to 6.5 change the unit of observation. Instead of one point per family averaged over everything, each point is a single named election, identified by a two-letter country code and a two-digit year: JP46 is Japan 1946, NZ96 New Zealand 1996.

Every election appears twice — once positioned by the country-level model and once by the district-level model, distinguished by font. The distance between a case’s two labels is therefore a visual reading of how much the choice of aggregation matters for that system. casesOfInterest below is the shared list of hand-picked elections; each of the three figures filters it to one family.

# Define a set of hand-picked elections of interest
# Each entry is a country code plus a two-digit year. The same vector serves
# Figures 6.3, 6.4 and 6.5, which each filter it down to one system family.
casesOfInterest <- c(
  "AU49", "AR46", "AR15",
  "BE14", "BO06", "BR14", "CL13", "EC98", "EC02",
  "FR12", "SL56", "GR52", "GR15", "IE97",
  "JP46", "JP47", "JP47", "MX79", "MX12", "MN08",
  "NZ93", "NZ96", "ZA99", "VE93", "VE05",
  "VE10", "AT02", "SN78", "ZA89"
)

# Identify relevant variables to extract
predictorsOfInterest <- c(
  "avg.M", "new.nvotes", "pool_level",
  "ballot_type", "threshold", "formula",
  "family", "totalEff.hat", "pers.hat", "id2",
  "expanded.family"
)

# Filter dataset to include only selected majoritarian elections
tempData <- realCases %>% dplyr::filter(family == "maj" & is.element(id2, casesOfInterest))

# Set up the I–I plot space
with(tempData, plot(
  x = totalEff.hat, y = pers.hat,
  xlab = "",
  ylab = "",
  xlim = c(0, 4), ylim = c(0.3, 0.6) #
  , axes = F, type = "n"
))

# Add axis titles
mtext("Interparty dimension", side = 1, line = 2)
mtext("(Total Duvergerian Effect)", side = 1, line = 3, cex = 0.8)
mtext("Intraparty dimension", side = 2, line = 3)
mtext("(Personalism Score)", side = 2, line = 2, cex = 0.8)
# Add axes
axis(1)
axis(2)
# Y-axis quadrant labels (vertical space of party orientation)
text(
  x = rep(2, 4), y = c(0.33, 0.41, 0.49, 0.57),
  labels = c(
    "Highly Centralized",
    "Party-Oriented",
    "Candidate-Oriented",
    "Highly Individualistic"
  ),
  col = "gray"
)
# X-axis quadrant labels (horizontal space of interparty permissiveness)
text(
  x = c(0.5, 1.25, 2.75, 3.5), y = rep(0.45, 4),
  labels = c(
    "Very Weak",
    "Permissive",
    "Constraining",
    "Very Strong"
  ),
  col = "gray", srt = 90
)
# Add COUNTRY-BASED GBM predictions using bold labels (font=4)
with(tempData, text(
  xy.coords(
    jitter(totalEff.hat, amount = 0.1),
    jitter(pers.hat, amount = 0.01)
  ),
  labels = id2, col = "black", cex = 0.8, font = 4
))
# Add DISTRICT-BASED GBM predictions using normal font
with(tempData, text(
  xy.coords(
    jitter(totalEff.hat.district, amount = 0.1),
    jitter(pers.hat.district, amount = 0.01)
  ),
  labels = id2, col = "black", cex = 0.8
))
# Add legend clarifying which label corresponds to which model
legend("bottomleft",
  legend = c("Country-based GBM", "District-based GBM"),
  text.font = c(4, 1), bty = "n"
)

For example, Argentina 1946 and Japan 1946 are positioned in the upper-left quadrant, characterized by relatively low Total Duvergerian Effect (TDE) scores and high Average Personalism (AP) scores, reflecting limited interparty constraints but strong candidate personalism. In contrast, systems like New Zealand 1996 and South Africa 1989, which later reformed away from SMDP, cluster toward the more constraining end of the interparty dimension. The plot also captures electoral reforms that shift incentives, such as Japan’s 1996 transition or Mongolia’s use of a two-round multiple non-transferable vote (TR-MNTV) system in 2008, illustrating how changes in electoral rules translate into meaningful movements within the I–I space.

Overall, Figure 6.3 demonstrates that majoritarian systems are not monolithic; their variation across countries and over time produces diverse interparty and intraparty incentives.

Figure 6.4: Select Proportional Electoral Systems

The following code generates Figure 6.4, which maps a selected set of proportional representation (PR) systems in the Interparty–Intraparty (I–I) space using predictions from both country-level and district-level GBM models. While PR systems generally produce low Total Duvergerian Effect (TDE) scores—indicating permissive incentives for party entry—this figure highlights significant variation in Average Personalism (AP) across PR systems.

# Subset dataset to proportional systems from the predefined cases of interest
# Two conditions joined by & : the system must be proportional AND its id must
# appear in the hand-picked list defined for Figure 6.3.
tempData <- realCases %>% dplyr::filter(family == "pr" & is.element(id2, casesOfInterest))
# Manually adjust IDs to differentiate between repeated entries for Greece 2015
# Greece changed from a closed to an open list within 2015, so the same country
# and year legitimately appears twice. The trailing [1] picks the FIRST of the
# two matching rows and relabels it, so the two can be told apart on the plot.
tempData$id2[tempData$id2 == "GR15"][1] <- "GR15a"
tempData$id2[tempData$id2 == "GR15"][1] <- "GR15b"

# Create the blank plot canvas
with(tempData, plot(
  x = totalEff.hat, y = pers.hat,
  xlab = "",
  ylab = "",
  xlim = c(-0.05, 4), ylim = c(0.3, 0.6),
  axes = F, type = "n"
))
# Add axis titles
mtext("Interparty dimension", side = 1, line = 2)
mtext("(Total Duvergerian Effect)", side = 1, line = 3, cex = 0.8)
mtext("Intraparty dimension", side = 2, line = 3)
mtext("(Personalism Score)", side = 2, line = 2, cex = 0.8)
# Draw axes
axis(1)
axis(2)
# Add quadrant labels to guide interpretation
text(
  x = rep(2, 4), y = c(0.33, 0.41, 0.49, 0.57),
  labels = c(
    "Highly Centralized",
    "Party-Oriented",
    "Candidate-Oriented",
    "Highly Individualistic"
  ),
  col = "gray"
)
text(
  x = c(0.5, 1.25, 2.75, 3.5), y = rep(0.45, 4),
  labels = c(
    "Very Weak",
    "Permissive",
    "Constraining",
    "Very Strong"
  ),
  col = "gray", srt = 90
)
# Add COUNTRY-LEVEL predictions (bold font)
with(tempData, text(xy.coords(totalEff.hat, pers.hat),
  labels = id2, col = "black", cex = 1, font = 4
))
# Add DISTRICT-LEVEL predictions (jittered, regular font)
with(tempData, text(
  xy.coords(
    jitter(totalEff.hat.district, amount = 0.1),
    jitter(pers.hat.district, amount = 0.01)
  ),
  labels = id2, col = "black", cex = 0.8
))
# Add a legend to clarify font type
legend("bottomleft",
  legend = c("Country-based GBM", "District-based GBM"),
  text.font = c(4, 1), bty = "n"
)

For example, closed-list PR (CLPR) systems such as Argentina 2015 (AR15) and Senegal 1978 (SN78) are positioned near the lower end of the personalism axis, while open-list PR (OLPR) systems like Brazil 2014 (BR14) and flexible-list systems like Austria 2002 (AT02) show much higher personalism scores. The figure also includes Greece 2015 twice (GR15a and GR15b) to illustrate how a shift from a closed to an open list within the same year increased personalism.

These differences primarily stem from ballot type and district magnitude. Figure 6.4 demonstrates that although PR systems are permissive in terms of interparty competition, their design choices can create substantial variation in intraparty incentives, with ballot structure playing a crucial role.

How to read Figures 6.3 through 6.5. All three share a frame, so they can be read as one argument in three parts. Within each, look for vertical spread: elections of the same family sitting at very different heights mean the family label does not determine intraparty incentives. Greece 2015, appearing twice in Figure 6.4 because the country switched list type mid-year, is the cleanest illustration — same country, same year, same magnitudes, different height.

Then compare across the three figures. Majoritarian cases should occupy the constraining end of the horizontal axis, PR cases the permissive end, and mixed systems should scatter between them. Finally, for any given election, note how far apart its two labels sit: a large gap means district-level and country-level models disagree about that system, which is exactly the methodological point made about mixed systems in Figure 6.5.

Figure 6.5: Select Mixed-Member Electoral Systems

Figure 6.5 visualizes a selection of mixed-member electoral systems in the Interparty–Intraparty (I–I) space, overlaying predictions from both country-level and district-level GBM models. This figure highlights the wide variation mixed systems exhibit depending on how their proportional and majoritarian components are combined.

# Filter realCases to include only mixed electoral systems from the selected cases
tempData <- realCases %>% dplyr::filter(family == "mixed" & is.element(id2, casesOfInterest))

# Initialize the Interparty–Intraparty (I–I) plot canvas
with(tempData, plot(
  x = totalEff.hat, y = pers.hat,
  xlab = "",
  ylab = "",
  xlim = c(0, 4), ylim = c(0.3, 0.6),
  axes = F, type = "n"
)) # Create blank plotting area

# Add axis and quadrant labels, labels the x-axis and y-axis and draws the coordinate axes.
mtext("Interparty dimension", side = 1, line = 2)
mtext("(Total Duvergerian Effect)", side = 1, line = 3, cex = 0.8)
mtext("Intraparty dimension", side = 2, line = 3)
mtext("(Personalism Score)", side = 2, line = 2, cex = 0.8)
# Draw axis ticks
axis(1)
axis(2)
# Add qualitative quadrant labels on the y-axis (personalism dimension)
text(
  x = rep(2, 4), y = c(0.33, 0.41, 0.49, 0.57),
  labels = c(
    "Highly Centralized",
    "Party-Oriented",
    "Candidate-Oriented",
    "Highly Individualistic"
  ),
  col = "gray"
)
# Add qualitative quadrant labels on the x-axis (interparty constraint dimension)
text(
  x = c(0.5, 1.25, 2.75, 3.5), y = rep(0.45, 4),
  labels = c(
    "Very Weak",
    "Permissive",
    "Constraining",
    "Very Strong"
  ),
  col = "gray", srt = 90
)
# Plot country-based GBM predictions (bold font, slight jitter to avoid overlapping text)
with(tempData, text(
  xy.coords(
    jitter(totalEff.hat, amount = 0.055),
    jitter(pers.hat, amount = 0.015)
  ),
  labels = id2, col = "black", cex = 0.8, font = 4
))
# Add district-based GBM predictions with regular font and more jitter
with(tempData, text(
  xy.coords(
    jitter(totalEff.hat.district, amount = 0.1),
    jitter(pers.hat.district, amount = 0.01)
  ),
  labels = id2, col = "black", cex = 0.8
))
# Add legend explaining the two labeling styles
legend("bottomleft",
  legend = c("Country-based GBM", "District-based GBM"),
  text.font = c(4, 1), bty = "n"
)

Some mixed systems, like Mexico 1979 (MX79) and Venezuela 2010 (VE10), cluster toward the more constraining end of the interparty dimension, reflecting stronger limitations on party entry. Others, such as Venezuela 1993 (VE93) and New Zealand 1996 (NZ96), fall closer to the permissive side, indicating more open competition. These differences arise from variations in vote pooling rules and the relative weighting of electoral tiers across countries.

The figure also shows discrepancies between country-level and district-level model predictions. For instance, Venezuela 2005 (VE05) and Mexico 2012 (MX12) shift positions depending on whether national configurations or district-level components inform the estimates. This illustrates an important methodological insight: without fine-grained district data, assessments of mixed systems’ incentives may vary substantially.

Figure 6.2: Placing Democracies on the I–I Space

The following code produces Figure 6.2, which maps only democratic electoral systems in the Interparty–Intraparty (I–I) space using country-level GBM predictions. By excluding elections held under anocratic or hybrid regimes, the plot focuses on how institutional design shapes electoral incentives within consolidated democracies.

Out of sequence. Figure 6.2 is produced here, at the end of the file, after Figures 6.3, 6.4 and 6.5. It appears earlier in the printed chapter.

# Filter to democratic cases, then group by electoral system family
# Same construction as Figure 6.1b, with one extra filter. Comparing the two
# figures answers the question of whether the pattern in 6.1 was driven by
# authoritarian or hybrid regimes holding elections under permissive rules.
mean_Democracies <- realCases %>%
  dplyr::filter(democ == 1) %>% # Keep only cases coded as democratic (Polity2 > 6)
  group_by(expanded.family) %>% # Group data by electoral system family
  summarize(
    mean.pers.hat = mean(pers.hat, na.rm = T) # Compute mean personalism score
    , mean.totalEff.hat = mean(totalEff.hat, na.rm = T) # Compute mean Total Duvergerian Effect
    , sd.pers.hat = sd(pers.hat, na.rm = T) # Compute standard deviation of personalism
    , sd.totalEff.hat = sd(totalEff.hat, na.rm = T)
  ) # Compute standard deviation of TDE

# Drop systems without defined family classification
mean_Democracies <- mean_Democracies[mean_Democracies$expanded.family != "missing", ]

# Create lowercase acronyms for selected PR subtypes
mean_Democracies$expanded.family.acronyms <- mean_Democracies$expanded.family
mean_Democracies$expanded.family.acronyms <- car::recode(
  mean_Democracies$expanded.family.acronyms,
  "'CLPR'='clpr';
                                          'FLPR'='flpr';
                                          'fort.PR'='fort.pr';
                                          'FrLPR'='frlpr';
                                          'OLPR'='olpr';
                                          'STV'='stv'"
)


# Initialize I–I plot space
par(mar = c(4, 4, 1, 2))
# with (mean_realCases, plot (x=mean.lsq.hat, y=mean.pers.hat
with(mean_Democracies, plot(
  x = mean.totalEff.hat, y = mean.pers.hat,
  xlab = "",
  ylab = "",
  xlim = c(0, 2.5), ylim = c(0.3, 0.6),
  axes = F, type = "n"
))
# Add axis labels and quadrant descriptors
mtext("Interparty dimension", side = 1, line = 2)
mtext("(Total Duvergerian effect)", side = 1, line = 3, cex = 0.8)
# mtext ("(average weighted variance)", side=1, line=3, cex=0.8)
# mtext ("(average disproportionality score)", side=1, line=3, cex=0.8)
mtext("Intraparty dimension", side = 2, line = 3)
mtext("(Personalism score)", side = 2, line = 2, cex = 0.8)
axis(1)
axis(2)
# Places qualitative interpretive labels on the plot axes,
text(
  x = rep(1.25, 4), y = c(0.33, 0.41, 0.49, 0.57),
  labels = c("Highly Centralized", "Party-Oriented", "Candidate-Oriented", "Highly Individualistic"),
  col = "gray"
)
text(
  x = c(0.1, 0.75, 1.75, 2.4), y = rep(0.45, 4),
  labels = c("Very Weak", "Permissive", "Constraining", "Very Strong"),
  col = "gray", srt = 90
)
# with (mean_realCases, segments (x0=mean.lsq.hat
#                                 , x1=mean.lsq.hat
#                                 , y0=(mean.pers.hat*(-1)+1) + sd.pers.hat
#                                 , y1=(mean.pers.hat*(-1)+1) - sd.pers.hat
#                                 , lwd=3, col="grey"))
# with (mean_realCases, segments (x0=mean.lsq.hat + sd.lsq.hat
#                                 , x1=mean.lsq.hat - sd.lsq.hat
#                                 , y0=(mean.pers.hat*(-1)+1)
#                                 , y1=(mean.pers.hat*(-1)+1)
#                                 , lwd=3, col="grey"))
# Draw ellipses and family labels
for (i in 1:nrow(mean_Democracies)) {
  # if ( is.element (mean_Democracies$expanded.family[i], c("AV", "SNTV") ) ) { next } else {
  points(
    ellipse(
      x = 0,
      centre = c(mean_Democracies$mean.totalEff.hat[i], mean_Democracies$mean.pers.hat[i]),
      scale = c(mean_Democracies$sd.totalEff.hat[i] / 2, mean_Democracies$sd.pers.hat[i] / 2)
    ),
    col = "gray", lwd = 1, lty = 2, type = "l"
  )
  # }
}
# Places electoral family acronyms at their respective locations in the I–I space.
with(
  mean_Democracies,
  text(xy.coords(mean.totalEff.hat, mean.pers.hat),
    labels = expanded.family.acronyms, col = color.scheme.2
  )
)
# Adds a legend explaining the color coding used for different electoral system types.
legend("bottomright",
  text.col = c("black", "slategray", "black"),
  legend = c("MAJORITARIAN SYSTEMS", "MIXED SYSTEMS", "proportional systems"),
  bty = "n"
)

Figure 6.2 locates only plausibly democratic electoral systems on the Interparty–Intraparty (I–I) Space, using predictions from a GBM model trained on country-level simulations. By excluding electoral contests held in anocracies and hybrid regimes, the figure aims to avoid distortions that may result from irregular or manipulated electoral configurations. The placement of each system reflects the average Total Duvergerian Effect (TDE) and Average Personalism (AP) scores, based on actual component rules. Dashed ellipses capture one standard deviation around these averages. As the chapter notes, removing anocratic elections does not dramatically alter the overall structure of the I–I Space, but some systems—such as limited nomination (LN), limited vote (LV), and multiple non-transferable vote (MNTV)—shift slightly. These adjustments lend confidence to the general patterns observed and reaffirm the utility of the I–I Space in mapping the incentives embedded in democratic electoral rules.

Moving Forward

This chapter has demonstrated how electoral systems, through their complex combinations of component rules, shape incentives for interparty and intraparty politics. By simulating thousands of elections and applying advanced machine learning models, we were able to locate real-world electoral systems within the Interparty–Intraparty (I–I) space. Our findings highlight substantial variation in how systems constrain or enable party competition (Total Duvergerian Effect, TDE) and how they promote either party-centered or candidate-centered politics (Average Personalism, AP). The diversity observed among mixed-member systems further illustrates the importance of considering the interplay between tiers and the granularity of data available for analysis.

Moving forward, the next part of the book (Part III) shifts from measuring incentives to exploring their consequences. Chapters 7 through 9 examine how interparty incentives affect the effective number of parties, the ideological distribution of party platforms, and the congruence between citizens’ preferences and their representatives’ positions. These chapters use empirical data to assess the real-world political outcomes linked to the incentive structures mapped here. Understanding this connection between electoral rules and party system dynamics is critical for scholars and practitioners aiming to design systems that foster desired political behaviors.

Session Information

Everything above was produced by a single run of this file. The two tables below record the environment that run took place in: the version of R, the machine and operating system, and the version of every package this chapter loads. They are generated while the page is being built, so they always describe the page you are reading rather than some earlier run.

This matters more than it may appear. R packages change: default arguments get revised, estimation routines are rewritten, and a number computed under one version of a modelling package is not guaranteed to reappear under another. Recording the versions is what allows a reader who gets a different answer to tell that apart from using a different version.

The R environment that produced this page
Setting	Value
R version	R version 4.5.1 (2025-06-13)
Platform	aarch64-apple-darwin20
Operating system	macOS Tahoe 26.5.2
Collation	en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
Time zone	America/Chicago
Pandoc	3.10.1
Page built on	30 July 2026

Packages this chapter attaches, with the version used here
Package	Version
car	3.1-3
carData	3.0-5
caret	7.0-1
doMC	1.3.8
dplyr	1.1.4
ellipse	0.5.0
forcats	1.0.0
foreach	1.5.2
gbm	2.2.2
ggplot2	3.5.2
ggpubr	0.6.3
gridExtra	2.3
gtools	3.9.5
iterators	1.0.14
lattice	0.22-7
lubridate	1.9.4
magrittr	2.0.3
mgcv	1.9-3
mixtools	2.0.0.1
nlme	3.1-168
openxlsx	4.2.8.1
purrr	1.2.2
readr	2.1.5
rpart	4.1.24
stringr	1.5.1
tibble	3.3.0
tidyr	1.3.1
tidyverse	2.0.0

The environment used for the published book

This chapter did not render its session-information block when the published results were produced, so no record of that original environment survives. The tables above describe this run only.

Where a record does survive, it is reproduced in full in RECOVERED-SessionInformation.md, alongside these files.

06. Place Real World Electoral Systems in the I-I Space – validation