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 11, 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 11 leverages variation across systems in incentives for intraparty politics — as captured by Average Personalism (AP) — to look at the extent to which elected officials report providing constituency service. It reflects on the advantages and disadvantages of different indicators of constituency service, then focuses on the PARTIREP dataset to construct a set of legislators’ constituency service scores. Empirical tests show that as AP rises legislators are more likely to devote time to providing constituency service.

This chapter builds its own outcome variable, and that is most of the work. Elsewhere in the book the dependent variable arrives ready-made — a count of bills, an effective number of parties. Here there is no single survey question asking “how much constituency service do you provide?”. Instead PARTIREP asks legislators about a range of related activities, and the chapter has to combine those answers into one measure.

The route runs: inspect the candidate survey items and how they correlate → extract a common factor from them by ordered factor analysis → use the resulting factor score as the outcome. Only then does the familiar machinery appear, with AP predicting that score.

So the first half of this file is measurement, not analysis. Figures 11.1 and 11.2 document the measure itself; Table 11.2 and Figures 11.3 and 11.4 are the results.

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. Conventions used throughout:

<- is R’s assignment arrow; x <- 5 stores 5 under the name x.
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).
A model formula like y ~ x reads “y is explained by x”.

Setup

Nothing on this page is cached. knitr can store a slow block’s results and reuse them on the next knit (the cache=TRUE chunk option). That speeds up editing, but it invalidates a stored result only when that block’s own code changes — not when a function or dataset it depends on is edited elsewhere in the file. A stale cache therefore produces output that looks current and is not, and one such cache did hide a real correction to this project’s code. Every block here runs on every knit. The file takes longer to compile; what it prints is always what the code currently says.

# Load relevant libraries
library(readxl) # reads .xlsx spreadsheets
library(tidyverse) # bundle of data-handling packages (dplyr, ggplot2, ...)
library(labelled) # handles the value labels carried by Stata (.dta) files
library(car) # recode() and other utilities
library(kableExtra) # formats knitr tables for HTML output
library(psych) # factor analysis -- used to build the constituency service measure
library(MCMCpack) # Bayesian estimation; supplies the ordered factor model
library(gbm) # gradient boosting machines: supply the AP and TDE predictions
library(mice) # multiple imputation; pools results across the five datasets
library(fixest) # fast fixed-effects estimation; supplies feols() below
library(sjPlot) # model plots
library(sandwich) # robust variance estimators
library(mixtools) # supplies rmvnorm(), used to simulate Figure 11.3

# Function to extract estimates for table 7.4
extract_Estimates <- function(model_name) {
  out <- ifelse(model_name$p.value < 0.01, paste0(round(model_name$estimate, 3), "***", " (", round(model_name$std.error, 3), ")"),
    ifelse(model_name$p.value < 0.05, paste0(round(model_name$estimate, 3), "**", " (", round(model_name$std.error, 3), ")"),
      ifelse(model_name$p.value < 0.1, paste0(round(model_name$estimate, 3), "*", " (", round(model_name$std.error, 3), ")"),
        paste0(round(model_name$estimate, 3), " (", round(model_name$std.error, 3), ")")
      )
    )
  )
  names(out) <- model_name$term
  out <- data.frame(out)
  out$variable <- row.names(out)
  return(out)
}

rubin_out <- function(model.list) {
  # Calculate Beta Bar (the coefficients)
  betabar <- Reduce("+", lapply(model.list, coef)) / length(model.list)
  # Calculate Uppsi. It is used in the calculations of SE's below.
  uppsi <- Reduce("+", lapply(model.list, function(x) (coef(x) - betabar) %*% t(coef(x) - betabar))) / (length(model.list) - 1)
  omegacl <- Reduce("+", lapply(model.list, function(x) vcov(x))) / length(model.list)
  # Calculate the final variance
  varR <- omegacl + as.numeric(uppsi) + as.numeric(uppsi) / length(model.list)
  # Return the results
  return(list(betas = betabar, vcov = varR))
}

Introduction

Chapter 11 explores the link between electoral system incentives and the extent to which legislators engage in constituency service—the activities legislators undertake to assist constituents with non-policy issues, such as navigating bureaucracy or securing local benefits. These behaviors are key mechanisms for building personal votes, especially in systems where intraparty competition encourages candidates to distinguish themselves.

Using survey data from the Parliamentary Elites of Latin America (PELA) and PARTIREP projects, the chapter constructs a latent constituency service score through ordinal factor analysis. This score synthesizes various self-reported measures of legislators’ constituency activities, providing a comprehensive measure of their engagement beyond legislative policymaking.

The core empirical analysis connects this constituency service score to Average Personalism (AP), our measure of electoral incentives for personal vote-seeking. Results indicate a positive and meaningful association: legislators operating in systems with higher AP scores tend to perform constituency service more frequently, even after controlling for individual, party, and country-level factors. The findings also highlight differences between national and regional legislators, with national MPs showing stronger responsiveness to AP incentives.

This file contains calls to datasets used in Chapter 11 (Constituency Service), as well as the code necessary to produce all graphs in the chapter.

Data processing

Data from PARTIREP project

As mentioned in Chapter 11, we extensively use data from the PARTIREP project, which surveys members of European national and subnational parliaments. The PARTIREP researchers have worked in 15 countries and 73 national and local assemblies. The following snippet loads the data (you can find a link to these data in our website) and produces a list of parliaments and survey years (each country+two digit code corresponds to an actual parliament; 00 are national parliaments, other numbers correspond to regional parliaments; the year under each country+two digit code conveys the year in which data were gathered.) At the end of the following snippet we use the sjlabelled package to extract the “label” attribute from Partirep’s Parliament indicator. This attribute makes observations readily recognizable as belonging to specific parliaments.

Along with PARTIREP, we make extensive use of a list of district magnitudes of regional parliaments that we collected online from various Wiki pages and other sources. These data (Reg.Parls in the snippet below) are loaded and merged with PARTIREP.

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/ch11/... 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.

# Set the path to Chapter 11 materials

# the short file names below. This is the one line to edit for your machine.


# Load district-level GBM objects
# Fitted models from Chapter 6. Nothing is estimated here -- they are only
# asked to predict AP and TDE for the systems the PARTIREP legislators sit in.
# The folder is spelled "RData/", matching every other chapter. It used to be
# "Rdata/" here; see the caveat above on why that mattered.
load(file = "data/shared/AP_district_objects_t515_d7.RData")
load(file = "data/shared/TDE_district_objects_t305_d15.RData")

# $optimalGBM holds the best-performing fitted model from each of five training
# runs, so each of these is a list of five models.
optimalGBMIntra <- totalAP.objects$optimalGBM
optimalGBMInter <- totalEffENP.objects$optimalGBM

# Load real systems scores
# ---------------------------------------------------------------------------
# This chapter previously read RealSystems_Scores_June_15_2023.RData, a
# 2023-era file, while chapters 7 to 10 read the later Aug 2024 file. The split
# was an oversight rather than a decision, and switching was verified to be
# costless here for a specific reason: this chapter does NOT use the stored TDE
# and AP scores at all. It recomputes them from the GBM models, using the
# electoral-rule variables carried by PARTIREP itself (the ES_* columns). The
# only thing it takes from this file is district magnitude, M -- and M is
# identical across the two files for every country-year the chapter keeps
# (62 shared rows, 0 differences), with no change in coverage.
#
# REQUIRES: a copy of RealSystems_Scores_GBM_Aug_2024.RData in this chapter's
# RData folder. It was never distributed here, so copy it from, for example,
# 07_The_Size_of_the_Party_System/Datasets/RData/ .
# ---------------------------------------------------------------------------
load("data/shared/RealSystems_Scores_GBM_Aug_2024.RData")

# Partirep
# The PARTIREP legislator survey, in Stata format. haven::read_dta() preserves
# Stata's value labels, which is why the `labelled` package is loaded above --
# the survey items arrive as labelled numbers rather than plain text.
Partirep <- haven::read_dta(file = "data/ch11/dta/PartiRep.dta")

# Regional parliaments M
# District magnitudes for regional assemblies, needed because PARTIREP covers
# both national and regional legislators. read_xlsx() returns a tibble; the
# [, 1:4] keeps the first four columns and as.data.frame() converts the result
# to an ordinary data frame.
Reg.Parls <- as.data.frame(read_xlsx("data/ch11/xlsx/Assemblies.xlsx", sheet = 1)[, 1:4])

Outcome variables

In Chapter 11, we consider nine items that come from question 15 in PARTIREP as indicators of “constituency service”. In PARTIREP, these items are preceded by the following text: “Mentioned below are some of the many different things that Members of Parliament do to keep in touch with constituents. For each one, do you actually do it outside election campaign periods, might you do it, or would you never?”. The nine items are:

Indicator	Description
`V015_1`	attending (or sending out letters on the occasion of) weddings, wedding anniversaries, and funerals in your local area
`V015_2`	meeting with (small parties of) constituents in their private home to talk about their wants and needs
`V015_3`	giving lectures and speaking at debate nights
`V015_4`	sending out a personal newsletter and direct mailing
`V015_5`	holding surgeries
`V015_6`	advertising your constituency work services (e.g. in newspaper ads or by visiting neighbourhoods)
`V015_7`	publicizing your successes in attracting business and obtaining government grants for the local area
`V015_8`	meeting local businesses and action groups featuring in the local media
`V015_9`	featuring in the local media

MPs are asked to provide one of six possible answers to such questions: “Actually do it at least once a week” (1), “Actually do it at least once a fortnight” (2), “Actually do it at least once a month” (3), “Actually do it at least every three months” (4), “might do it” (5), “would never” (6). In other words, higher values correspond to less constituency service. PARTIREP uses two codes for missing values: 99999 and 88888.

# Find NA codes and turn into NA
Partirep <- Partirep %>% dplyr::mutate(across(where(is.numeric), ~ dplyr::na_if(., 99999)))
Partirep <- Partirep %>% dplyr::mutate(across(where(is.numeric), ~ dplyr::na_if(., 88888)))

# Recode tier
Partirep$Tier <- car::recode(Partirep$Tier, "4=0")
Partirep$Parliament.tier <- paste0(Partirep$Parliament.ID, "-0", Partirep$Tier)

# Inspect potential outcome variables separately
outcomeVariables <- c(
  "V015_1", "V015_2", "V015_3",
  "V015_5", "V015_8", "V015_4", "V015_6",
  "V015_7", "V015_9"
)
# moreOutcomeVariables <- c("V003_1","V003_2","V002_3","V016","V019","V044","V055")
moreOutcomeVariables <- c("V044")

# Select outcomes
Outcomes <- dplyr::select(Partirep, which(names(Partirep) %in% outcomeVariables))

# Values 7 and 8 should be NA in V015
for (i in 1:length(outcomeVariables)) {
  val_label(Outcomes[, grep(outcomeVariables[i], colnames(Outcomes))], 7) <- NULL
  val_label(Outcomes[, grep(outcomeVariables[i], colnames(Outcomes))], 8) <- NULL
  na_values(Outcomes[, grep(outcomeVariables[i], colnames(Outcomes))]) <- c(7, 8)
  Outcomes[, grep(outcomeVariables[i], colnames(Outcomes))] <- user_na_to_na(Outcomes[, grep(outcomeVariables[i], colnames(Outcomes))])
}

# Select variables
Outcomes <- Outcomes %>% dplyr::select(V015_1, V015_2, V015_3, V015_5, V015_8, V015_4, V015_6, V015_7, V015_9)

# Covert to numeric
Outcomes <- lapply(Outcomes, as.numeric)

# Create a dataset with otucomes
Outcomes <- as.data.frame(do.call(cbind, Outcomes))

# We need to change the order of the labels, so that higher numbers (6) correspond to
# more frequency of constituency service
Outcomes <- as.data.frame(apply(Outcomes, 2, car::recode, "6=1; 5=2; 4=3; 3=4; 2=5; 1=6"))
Ambition <- dplyr::select(Partirep, which(names(Partirep) %in% moreOutcomeVariables))

# ambition gets 1 if the MP intends to stand in elections in next cycle
# 0 if not, or has not arrived at a decision
# NA otherwise
ambition <- car::recode(as.numeric(Ambition$V044), "3=0; 2=0; 4=NA; 5=NA")

# Recode Partirep variables that will be used later
Partirep$Rookie <- car::recode(Partirep$Rookie, "1=0; 2=1; 3=NA")
Partirep$Leadpos <- car::recode(Partirep$Leadpos, "1=0; 2=1; 3=NA")
Partirep$Sex <- car::recode(Partirep$Sex, "2=0")
Partirep$GOV <- car::recode(Partirep$GOV, "1=0; 2=1")
Partirep$REG <- car::recode(Partirep$REG, "1=0; 2=1")
Partirep$ambition <- ambition

Distribution of outcome variables

The following table simply summarizes each of the nine variables that we see as manifestations of constituency service:

outcomeTable <- Outcomes %>%
  mutate(Country = group_indices(., Partirep$Country)) %>%
  group_by(Country) %>%
  summarise(
    weddings = mean(V015_1, na.rm = T),
    meetings = mean(V015_2, na.rm = T),
    lectures = mean(V015_3, na.rm = T),
    surgery = mean(V015_5, na.rm = T),
    meetgreet = mean(V015_8, na.rm = T),
    newslttr = mean(V015_4, na.rm = T),
    advertise = mean(V015_6, na.rm = T),
    publicize = mean(V015_7, na.rm = T),
    localmedia = mean(V015_9, na.rm = T)
  )
outcomeTable$Country <- names(val_labels(Partirep$Country))
kable(outcomeTable, digits = 2) %>%
  kable_styling(font_size = 10, "striped")

Country	weddings	meetings	lectures	surgery	meetgreet	newslttr	advertise	publicize	localmedia
AUT	3.44	3.24	4.17	4.17	3.54	4.17	3.44	2.66	4.86
BEL	3.28	2.83	3.90	3.77	3.26	3.45	2.23	2.72	4.15
FRA	2.49	2.32	3.48	4.35	3.61	3.38	3.66	2.95	4.19
GER	3.13	2.85	4.57	4.11	3.92	4.26	3.31	3.03	5.06
HUN	2.24	3.55	4.12	4.16	3.31	2.83	4.12	3.00	4.47
IRE	4.62	3.56	2.97	5.18	4.38	3.15	4.09	3.58	5.24
ISR	3.33	3.72	4.26	3.39	1.92	2.92	2.82	1.97	3.37
ITA	2.95	3.15	4.71	4.81	3.81	3.72	3.23	2.49	4.61
NET	2.07	2.97	4.75	1.82	3.72	3.10	1.38	2.86	3.84
NOR	1.93	3.86	4.30	2.13	4.07	2.72	NaN	3.02	4.98
POL	2.10	2.27	3.84	5.64	3.59	2.53	4.00	2.96	4.27
POR	1.92	2.70	3.74	4.14	3.42	4.10	3.64	2.36	4.45
SPA	2.35	2.68	4.52	3.88	3.68	4.12	2.91	2.79	4.75
SWI	1.85	2.13	3.11	1.80	2.46	2.13	2.15	1.79	3.34
UNK	2.20	3.15	3.27	5.01	4.31	3.24	4.00	3.43	5.34

Most of these items refer to activities that require the MP’s physical presence among constituents, very likely involving trips to the MP’s actual circumscription (events, meetings, lectures, surgeries, meetgreet). The other items (newsletter, advertise, publicize, localmedia) should be more appropriately described as credit-claiming behavior, but they still concern advertising for constituency-oriented activities. In the pooled dataset (i.e., looking at the distribution of outcomes without breaking down the data at the parliament level), these nine items are indeed correlated. The following table reports Kendall’s tau correlation coefficients, a nonparametric measure well-suited for ordinal data, for all nine items.

# Compute the correlation matrix using Kendall's tau (appropriate for ordinal data)
corTable <- cor(Outcomes, use = "p", method = "kendall")

# Assign more readable row and column names to the correlation matrix
rownames(corTable) <- colnames(corTable) <- c(
  "events", # Attend/send letters for local events
  "meetings", # Meet with constituents at home
  "lectures", # Give lectures or talks
  "surgery", # Hold one-on-one problem-solving meetings
  "meetgreet", # Meet local businesses and action groups
  "newslttr", # Send personal newsletter or direct mail
  "advertise", # Advertise constituency work
  "publicize", # Publicize achievements
  "localmedia"
) # Appear in local media outlets


# Display the correlation matrix using kable with a descriptive caption and formatting
kable(corTable, caption = "Outcome correlations", digits = 2) %>%
  kable_styling(font_size = 12, "striped")

Outcome correlations
	events	meetings	lectures	surgery	meetgreet	newslttr	advertise	publicize	localmedia
events	1.00	0.32	0.16	0.24	0.24	0.18	0.20	0.28	0.20
meetings	0.32	1.00	0.27	0.23	0.29	0.23	0.24	0.25	0.24
lectures	0.16	0.27	1.00	0.20	0.35	0.36	0.17	0.19	0.37
surgery	0.24	0.23	0.20	1.00	0.34	0.32	0.41	0.32	0.34
meetgreet	0.24	0.29	0.35	0.34	1.00	0.30	0.30	0.45	0.47
newslttr	0.18	0.23	0.36	0.32	0.30	1.00	0.26	0.25	0.32
advertise	0.20	0.24	0.17	0.41	0.30	0.26	1.00	0.40	0.32
publicize	0.28	0.25	0.19	0.32	0.45	0.25	0.40	1.00	0.35
localmedia	0.20	0.24	0.37	0.34	0.47	0.32	0.32	0.35	1.00

And the following snippet produces Plot (a) in Figure 11.1, i.e., a “shaded box” version of Kendall’s tau correlation coefficients shown in the table immediately above.

# Set plot margins to allow room for axis labels
par(mar = c(4, 4, 0, 0))

# Generate a lower-triangle correlation matrix plot
cor.plot(corTable,
  numbers = T, upper = FALSE,
  main = "" # Kendall's tau correlation among constituency service items
  , show.legend = FALSE
)

# dev.print(pdf        # copies the plot to a the PDF file
# , width=9, height=7
# , "corTable-constituency-service.pdf")

Figure 11.1(a) confirms that all pairs of constituency service activities are positively correlated, supporting the idea that they capture a common latent dimension of MP behavior. Notably, meetgreet, localmedia, and publicize exhibit some of the strongest correlations, suggesting that personal visibility and interaction are central components of constituency service. Conversely, items like lectures and events show weaker links to other practices, indicating they may reflect more ceremonial or symbolic outreach. The overall structure justifies modeling these items as expressions of a single latent trait in the following ordinal factor analysis.

Ordered factor analysis

Rather than looking at each of the nine manifest indicators of constituency service one by one, in Chapter 11 we build a factor analytic model that allows us to infer where each MP falls along a continuous “constituency service” scale. Factor analysis is part of a family of statistical models known as “scaling” or “decomposition” techniques. They are very useful in “reducing the dimensionality” of a problem; you can see here that, rather than presenting nine different models, one per outcome, we can present a single model based on “constituency service” scores as the outcome of interest.

More importantly, we can think of “scaling techniques” as models that allow us to capture a latent, unobservable phenomenon. In this case, we are obviously interested in an MP’s propensity to engage in constituency service, but this theoretical concept is not easily measurable. What we can measure instead are many manifest indicators that suggest that constituency service is happening. The nine outcome variables that we have inspected are the manifest indicators, and the scaling techniques that we use allow us to “score” the extent of constituency service in which MPs engage.

We opted for a more complex Bayesian factor analysis model that takes into account the fact that manifest variables are in fact ordered categories. In other words, the numbers 1 through 6 with which we code responses to manifest indicators cannot be interpreted as capturing cardinal information. An MP that chooses 4, instead of 2, cannot be seen as “doubling” her constituency service effort. A model that appropriately accounts for ordinal data is thus required, and a Bayesian inferential framework makes it easier to estimate such a model. A good primer into scaling techniques is Johnson and Albert 1999. (Please note that running the following step takes a few minutes, even on a fast computer):

What the factor model does, and why it is ordered. Nine survey items each capture a piece of constituency service, and no single one is the concept. Factor analysis assumes an unobserved trait — call it “propensity to serve the constituency” — that drives a legislator’s answers to all nine, and recovers a score for each legislator on that hidden trait.

It has to be an ordered (ordinal) factor model because the answer scale is ranked but not evenly spaced. “At least once a week” through “would never” runs in a clear order, yet the gap between weekly and fortnightly is not the same quantity as the gap between “might do it” and “would never”. Treating those codes as numbers would assume a spacing the survey never established.

Note also the direction of the scale, stated in the text above: higher raw values mean less constituency service. Keep that in mind when reading the loadings and the figures.

This is Bayesian estimation via MCMC, so it explores the space of plausible answers by simulation rather than solving for one solution. seed = 1971 fixes those draws, so repeated runs give identical answers. The step is slow, and it is recomputed on every knit rather than stored between runs.

# Run Bayesian ordinal factor analysis (Quinn 2004), assuming 1 latent factor
# The model is identified under rotational invariance by forcing the factor loading for item V015_1 to be positive
# Why a constraint is needed: a factor model cannot by itself tell "high score
# means more service, all loadings positive" from the mirror image with every
# sign flipped -- both fit identically. Pinning one item's loading to be
# positive picks an orientation so the scores are interpretable.

ordFactAnal <- MCMCordfactanal(as.matrix(Outcomes),
  factors = 1,
  lambda.constraints = list("V015_1" = list(2, "+")),
  burnin = 1000 # Burn-in period for MCMC
  , mcmc = 10000 # Number of MCMC draws
  , thin = 100 # Thinning to reduce autocorrelation
  , verbose = FALSE,
  seed = 1971 # Reproducibility
  , lambda.start = 2 # Starting value for loading
  , store.scores = T
) # Save factor scores

## 
## 
## Acceptance rates:
##  V015_1 V015_2 V015_3 V015_5 V015_8 V015_4 V015_6 V015_7 V015_9
##    0.72   0.78   0.77   0.69   0.81   0.72   0.76    0.8   0.79

# On those settings: MCMC wanders toward the right answer, so the first 1,000
# draws (the burn-in) are discarded as unrepresentative. Of the 10,000 kept,
# thin=100 retains every hundredth, since consecutive draws resemble one
# another and add little information.

# Save MP-specific factor scores; these are the constituency service scores that we analyze in Chapter 11.
# The output stores one column per legislator (named "phi..."). grep() finds
# those columns and colMeans() averages each legislator's draws into a single
# score -- the outcome variable for the rest of the chapter.
ordFA.scores <- colMeans(ordFactAnal[, grep("phi", colnames(ordFactAnal))])

# Extract item loadings from the posterior draws
# A loading says how strongly one survey item reflects the underlying trait.
# Items with large loadings are good indicators of constituency service; items
# near zero contribute little.

ordFA.loadings <- ordFactAnal[, grep("Lambda", colnames(ordFactAnal))]
# The model stores two parameters per item -- a difficulty and a loading -- as
# columns named Lambda1.1, Lambda1.2, and so on. The '\\.2' pattern keeps only
# the ".2" columns, which are the loadings. The double backslash escapes the
# dot, which would otherwise match any character in a regular expression.
ordFA.loadings <- ordFA.loadings[, grep("\\.2", colnames(ordFA.loadings))]

# Rename for interpretability
colnames(ordFA.loadings) <- c(
  "weddings", "meetings", "lectures",
  "surgery", "meetgreet", "newslttr",
  "advertise", "publicize", "localmedia"
)

How to read Figure 11.2. The distribution of the newly built factor scores, broken out by country. Still measurement rather than result: it shows the outcome variable behaving sensibly before it is asked to do any work.

Look for spread both within countries — legislators in the same system differ, which is what makes an individual-level analysis meaningful — and between them, since systematic country differences are what AP will later be asked to explain.

Figure 11.2: Constituency Service Factor Scores

We produce a boxplot that breaks down scores by Parliament-tier.

# Boxplot of constituency service, with tier-Parliament as grouping category
Partirep$Tier <- car::recode(Partirep$Tier, "4=0") # Recode Tier 4 (list tier) as 0

# Create a unique parliament-tier identifier
Partirep$Parliament.tier <- paste0(Partirep$Parliament.ID, "-0", Partirep$Tier)

# Combine scores and tiers into a single data frame
boxplotDataTier <- data.frame(constservice = ordFA.scores, parliament = Partirep$Parliament.tier)

# Order the boxplot by the median constituency score in each parliament-tier
new_order_Tier <- with(boxplotDataTier, reorder(parliament, constservice, median, na.rm = T))
# Assign colors: white = SMD tier (Tier 1), gray = list tier (Tier 0)
boxplotColorTier <- rep("gray", length(levels(new_order_Tier)))
boxplotColorTier[grep("-01", levels(new_order_Tier))] <- "white"
# Highlight Germany’s Bundestag with black labels
col.labels <- rep("gray", length(levels(new_order_Tier)))
col.labels[grep("GER00", levels(new_order_Tier))] <- "black"

mean(boxplotDataTier$constservice[boxplotDataTier$parliament == "GER00-01"]) # avg constituency score in national parliament

## [1] 0.688

mean(boxplotDataTier$constservice[boxplotDataTier$parliament == "GER00-02"]) # avg constituency score in regional parliament

## [1] 0.254

# Helper to determine even/odd for alternating axis labels
is.odd <- function(x) x %% 2 == 1

# Plot configuration
par(mar = c(3, 4, 3, 1))
boxplot(boxplotDataTier$constservice ~ new_order_Tier, ylab = "Constituency service scores", xlab = "", axes = F, col = boxplotColorTier)
# Add axes and labels with alternating position and rotation
axis(2)
axis(1, at = c(1:length(levels(new_order_Tier)))[is.odd(c(1:length(levels(new_order_Tier))))], labels = FALSE)
axis(3, at = c(1:length(levels(new_order_Tier)))[!is.odd(c(1:length(levels(new_order_Tier))))], labels = FALSE)
# Add tier labels to avoid overlap
text(seq(1, length(levels(new_order_Tier)), by = 2), par("usr")[3] - 0.2,
  labels = levels(new_order_Tier)[is.odd(c(1:length(levels(new_order_Tier))))],
  srt = 45, pos = 1, xpd = TRUE, cex = 0.5,
  col = col.labels[is.odd(c(1:length(levels(new_order_Tier))))]
)
text(seq(2, length(levels(new_order_Tier)), by = 2), par("usr")[4] + 0.5,
  labels = levels(new_order_Tier)[!is.odd(c(1:length(levels(new_order_Tier))))],
  srt = 45, pos = 1, xpd = TRUE, cex = 0.5,
  col = col.labels[!is.odd(c(1:length(levels(new_order_Tier))))]
)
# Add legend distinguishing tiers
legend("bottomright",
  legend = c("SMD tier", "List tier"),
  pt.bg = c("white", "gray"),
  pch = 21, bty = "n"
)

# dev.print(pdf        # copies the plot to a the PDF file
#          , width = 9, height = 7
#          , "boxplot-constituency-service-parliament-tier-level.pdf")

Figure 11.2 displays boxplots of the constituency service factor scores across different parliament-tier combinations. Two patterns are clear. First, MPs elected via SMD tiers (white boxes) tend to report higher levels of constituency service than those from list tiers (gray boxes), consistent with theories that link single-member systems to personal vote cultivation. Second, even within countries, substantial variation exists—most strikingly in Germany, where Bundestag MPs elected in nominal tiers (GER00-01) score far higher than their list-tier counterparts (GER00-02). This figure supports the claim that electoral rules shape incentives for constituency work more strongly than the mere level of government.

Building AP scores from PARTIREP predictors

PARTIREP includes a number of macro variables (i.e., variables at the parliament-tier or parliament-district levels, rather than at the individual level) that describe the electoral system under which individual MPs were elected.

Indicator	Description
`Tier`	1, 2, or 3 (88888 not applicable, 99999 missing)
`Candidate`	1=at nominal level, 2=at proportional level, 3=both levels (88888 not applicable, 99999 missing)
`ES_Descr`	a: PR, b: STV, c: TR, d: Plurality
`ES_Div`	a: LR-Hare, b: LR-Droop, c: LR-Imperiali, d: Imperiali, e: D’Hondt, f: Mod SL, g: SL, h: equal proportions, i: Danish, j: STV - Droop, k: other
`ES_MM_seat`	1= list is at least partially compensatory?
`ES_MM_vote`	1= at least two votes?
`ES_Thres`	electoral threshold (%)
`ES_Vote_P`	Voters can vote for a party list (1=yes)
`ES_Vote_C`	Votes can vote for a candidate (1=yes)
`ES_Pool`	Are votes for candidates pooled? a: across party list, b: below party level, c: no
`ES_Ballot`	List type: a: open, b: flexible but rather open, c: flexible but rather closed, d: closed

# Select relevant electoral rule variables from the PARTIREP dataset
elecVars <- c(
  "Tier", "Candidate", "ES_Descr", "ES_Div",
  "ES_MM_seat", "ES_MM_vote", "ES_Thres", "ES_Vote_P",
  "ES_Vote_C", "ES_Pool", "ES_Ballot", "ES_Vote_max", "ES_Vote_pan", "REG"
)
ElectoralPredictors <- dplyr::select(Partirep, which(names(Partirep) %in% elecVars))

# Recodifications of variables so that they coincide with names in "realCases" dataset
# pool_level
# translation here is: "candidate"="no pooling"; "party"="across the party list"; "party_list"="below the party level". This might not be right, check PARTYREP manual
pool_level <- base::as.factor(ElectoralPredictors$ES_Pool)
levels(pool_level) <- c("party", "party_list", "candidate")

# ballot_type (turn into "closed", "open", "flexible")
ballot_type <- base::as.factor(ElectoralPredictors$ES_Ballot)
levels(ballot_type) <- c(NA, "closed", "flexible", "flexible", "open")

# new.nvotes (turn into One, TotalCandidates, TotalSeats, LessSeats)
# vote max counts the number of preferrential votes, we need to recode it
ElectoralPredictors$ES_Vote_max_new <- (ElectoralPredictors$ES_Vote_max)
ElectoralPredictors$ES_Vote_max_new[ElectoralPredictors$ES_Vote_max > 1 & ElectoralPredictors$ES_Descr == 3 & ElectoralPredictors$ES_Vote_pan == 1] <- 1
new.nvotes <- as.factor(ElectoralPredictors$ES_Vote_max_new)
levels(new.nvotes) <- c("One", NA, "TotalSeats", "TotalCandidates")

# Formula
# we need to get to:
# [1] plurality         abs majority      hare
# [4] dhondt            hagenbachbischoff <NA>
# [7] saintelague       modsaintlague     imperiali
# [10] droop

# ES_Descr defines majoritarian vs PR logic
# Map to basic family categories (first step in formula creation)
formula1 <- base::as.factor(ElectoralPredictors$ES_Descr)
levels(formula1) <- c("abs majority", "plurality", "PR", "PR")
formula1 <- as.character(formula1)

# Recode ES_Div (divisor formulas used in PR)
formula2 <- base::as.factor(ElectoralPredictors$ES_Div)
# NOTE: ES_Div includes "other", but does not include "hagenbachbischoff";
# I codify "other" as "hagenbachbischoff"
levels(formula2) <- c(
  NA, "dhondt", "droop",
  "hare", "modsaintlague",
  "hagenbachbischoff",
  "saintelague"
)
formula2 <- as.character(formula2)

# Final formula variable: for PR systems use the PR-specific rule, otherwise use plurality/majority type
formula <- ifelse(formula1 == "plurality" | formula1 == "abs majority", formula1, formula2)

# Use Candidate variable for mixed-member systems
# family is one of "maj", "mixed", "pr"
family1 <- base::as.factor(ElectoralPredictors$Candidate)
levels(family1) <- c("mixed", "mixed", "mixed", NA)
family1 <- as.character(family1)

# Use ES_Descr to define whether system is majoritarian or PR
family2 <- base::as.factor(ElectoralPredictors$ES_Descr)
levels(family2) <- c("maj", "maj", "pr", "pr")
family2 <- as.character(family2)

# Combine both sources of information
family <- ifelse(is.na(family1), family2, family1)

# Threshold (number between 0 and 17)
threshold <- ElectoralPredictors$ES_Thres
threshold <- ifelse(is.na(threshold), 0, threshold)

# Add all rules to ElectoralPredictors
ElectoralPredictors$pool_level <- pool_level
ElectoralPredictors$ballot_type <- ballot_type # has missing values
ElectoralPredictors$new.nvotes <- new.nvotes # has missing values
ElectoralPredictors$formula <- formula
ElectoralPredictors$family <- family
ElectoralPredictors$threshold <- threshold

PARTIREP does not have information on district magnitude. For legislators that compete in plurality-majoritarian systems, and for legislators that compete in single-member tier 1 districts in mixed-member electoral systems, we can easily add M=1 as a predictor. For other legislators, we will employ the average district magnitude that we obtain from our own dataset.

# List of countries covered in the PARTIREP data
countries2keep <- c(
  "Austria", "Belgium", "France", "Germany", "Hungary",
  "Ireland", "Israel", "Italy", "Netherlands", "Norway",
  "Poland", "Portugal", "Spain", "Switzerland", "United Kingdom"
)

# Load the realCases data and copy M into avg.M
realCases <- tibble::as_tibble(data2export)
realCases$avg.M <- realCases$M # We need avg.M below

# Keep post-2000 elections in countries of interest, add ISO-style country ID
partirepCases <- realCases %>%
  dplyr::filter(country %in% countries2keep & year > 2000) %>%
  dplyr::mutate(id = dplyr::recode(country,
    "Austria" = "AUT",
    "Belgium" = "BEL",
    "France" = "FRA",
    "Germany" = "GER",
    "Hungary" = "HUN",
    "Ireland" = "IRE",
    "Israel" = "ISR",
    "Italy" = "ITA",
    "Netherlands" = "NET",
    "Norway" = "NOR",
    "Poland" = "POL",
    "Portugal" = "POR",
    "Spain" = "SPA",
    "Switzerland" = "SWI",
    "United Kingdom" = "UNK"
  ))

# Create a unique identifier for each case: country-year
partirepCases$id_year <- paste(partirepCases$id, partirepCases$year, sep = "-")

Create average district magnitude for systems in dataset

# Create a country ID vector for the PARTIREP data
id <- base::as.factor(Partirep$Country)
levels(id) <- c("AUT", "BEL", "FRA", "GER", "HUN", "IRE", "ISR", "ITA", "NET", "NOR", "POL", "POR", "SPA", "SWI", "UNK")
id <- as.character(id)

# Generate identifiers for joining with realCases
ElectoralPredictors$id_year <- paste(id, Partirep$Start_term, sep = "-")
ElectoralPredictors$id <- id
ElectoralPredictors$year <- Partirep$Start_term

# In addition, ballot_type and new.nvotes have many missing values
# might be best to also get them from partirepCases
avg.magnitude <- ballot <- nvotes <- c()

# Loop through unique id-year combinations
for (i in 1:length(unique(ElectoralPredictors$id_year))) {
  country.year <- unique(ElectoralPredictors$id_year)[i]
  country <- ElectoralPredictors$id[ElectoralPredictors$id_year == country.year]
  year <- ElectoralPredictors$year[ElectoralPredictors$id_year == country.year]
  unique.country <- unique(country)
  unique.year <- unique(year)
  tmp.data <- partirepCases[partirepCases$id == unique.country, ]

  # Match M, ballot, and nvotes either directly or from the most recent past year
  if (is.element(country.year, tmp.data$id_year)) {
    tmp.avg.M <- tmp.data$avg.M[tmp.data$id_year == country.year]
    tmp.ballot <- as.factor(tmp.data$ballot_type[tmp.data$id_year == country.year])
    tmp.nvotes <- as.factor(tmp.data$new.nvotes[tmp.data$id_year == country.year])
  } else {
    tmp.data <- tmp.data[tmp.data$year < unique.year, ]
    tmp.avg.M <- tmp.data$avg.M[which.min(unique.year - tmp.data$year)]
    tmp.ballot <- as.factor(tmp.data$ballot_type[which.min(unique.year - tmp.data$year)])
    tmp.nvotes <- as.factor(tmp.data$new.nvotes[which.min(unique.year - tmp.data$year)])
  }

  # Expand values to match number of MPs
  avg.magnitude <- c(avg.magnitude, rep(tmp.avg.M, length(country)))
  ballot <- c(as.factor(ballot), rep(tmp.ballot, length(country)))
  nvotes <- c(as.factor(nvotes), rep(tmp.nvotes, length(country)))
}

# Adding avg.M requires bringing in information from our dataset
# Assign M = 1 for majoritarian and SMD systems
avg.M <- ifelse(ElectoralPredictors$ES_Descr == 1 | ElectoralPredictors$ES_Descr == 2, 1, avg.magnitude)

# But now we will eschew the avg.M values from our dataset for
# regional parliaments in favor of the M values that BC obtained
# in the various regional parliament wiki pages detailed in
# Assemblies.xlsx
avg.M <- ifelse(Partirep$REG == 1, Partirep$m, avg.M)

# Now add average magnitude to ElectoralPredictors
ElectoralPredictors$M <- avg.M

# ballot type
ballot <- as.character(ballot)
ElectoralPredictors$ballot_type <- as.character(ElectoralPredictors$ballot_type)
ElectoralPredictors$ballot_type <- ifelse(is.na(ElectoralPredictors$ballot_type), ballot, ElectoralPredictors$ballot_type)
ElectoralPredictors$ballot_type <- as.factor(ElectoralPredictors$ballot_type)

# Votes
ElectoralPredictors$new.nvotes <- as.character(ElectoralPredictors$new.nvotes)
ElectoralPredictors$new.nvotes <- ifelse(is.na(ElectoralPredictors$new.nvotes), "One", ElectoralPredictors$new.nvotes) # if value is missing for new.nvotes, we assume it is one; this actually coincides with nvotes, which comes from our own dataset
ElectoralPredictors$new.nvotes <- as.factor(ElectoralPredictors$new.nvotes)

# formula
ElectoralPredictors$formula <- as.factor(ElectoralPredictors$formula)

We now have in the R object ElectoralPredictors enough information to find all relevant covariates that, combined with the GBM predictive model, allow us to arrive at predicted values of TDE and AP to characterize the electoral environment in which each of the MPs in PARTIREP is located. The following code snippet loads the predictive GBMs and produces two variables: totalEff.hat (TDE) and pers.hat (AP).

The original code was using the old GBM. Here I am using the first GBM in the list just to make sure that the code runs correctly and produces results similar to the ones in the book.

# Gather relevant variables into a data.frame
Model.Data.Base <- data.frame(
  constituency.service = ordFA.scores # as.numeric(poly_model$scores)
  , party = as.factor(Partirep$Party),
  country = as.factor(Partirep$Country),
  sex = Partirep$Sex # 0=Female, 1=Male
  , rookie = Partirep$Rookie # 1=newly elected
  , leader = Partirep$Leadpos # 1=speaker, committee chair, PPG leader
  , govParty = Partirep$GOV # 1=government party
  , regionalParl = Partirep$REG # 1=regional party
  , ideology = as.factor(Partirep$P_ideol),
  ambition = Partirep$ambition,
  parl.tier = Partirep$Parliament.tier,
  M = ElectoralPredictors$M,
  ballot_type = ElectoralPredictors$ballot_type,
  Parliament.ID = Partirep$Parliament.ID
)

# Predict TDE and AP using each of the 5 imputed GBM models
Model.Data <- list()
# Drop observations where M is missing (e.g., regional parliaments in France)
for (i in 1:5) {
  # Interparty
  Model.Data.Base$TDE <- predict(optimalGBMInter[[i]], ElectoralPredictors, n.trees = optimalGBMInter[[i]]$n.trees)
  # Intraparty
  Model.Data.Base$AP <- predict(optimalGBMIntra[[i]], ElectoralPredictors, n.trees = optimalGBMIntra[[i]]$n.trees)
  # Put in a list
  Model.Data[[i]] <- Model.Data.Base
}

# Remove regional parliaments from France (we don't have data for M)
# These systems uses a complex system in which the party that wins 50% of the votes in the first round
# gets 1/4 of the seats and the other seats are distributed proportionally (5% threshold).
# Otherwise, there is a second round
# Source: https://fr.wikipedia.org/wiki/%C3%89lections_r%C3%A9gionales_fran%C3%A7aises_de_2010
for (i in 1:5) Model.Data[[i]] <- subset(Model.Data[[i]], !is.na(M))

Four models, and what the last one is for. The outcome throughout is constituency.service, the factor score built above.

m1 AP alone
m2 adds TDE, checking AP is not standing in for interparty incentives
m3 adds legislator characteristics: sex, rookie (first term), leader (party office), govParty (governing party), regionalParl (sits in a regional rather than national assembly)
m4 drops the regional-parliament control and instead restricts the sample to national legislators only, via subset(..., regionalParl == 0)

m3 and m4 are two ways of handling the same worry. m3 keeps everyone and adjusts for tier; m4 sets the regional legislators aside entirely. Agreement between them means the finding does not depend on how that mixed sample is treated.

| country adds country fixed effects, so comparisons are made within a country. cluster = 'parl.tier' allows for correlated errors among legislators sitting in the same chamber, rather than assuming each responds independently.

# Initialize model lists
# Four containers, each of which will hold five fitted models -- one per
# imputation of the GBM predictions.
m1_l <- list()
m2_l <- list()
m3_l <- list()
m4_l <- list()


# Fit four models on each imputation
for (i in 1:5) {
  m1_l[[i]] <- feols(constituency.service ~ AP | country, cluster = "parl.tier", data = Model.Data[[i]])
  m2_l[[i]] <- feols(constituency.service ~ TDE + AP | country, cluster = "parl.tier", data = Model.Data[[i]])
  m3_l[[i]] <- feols(constituency.service ~ TDE + AP + sex + rookie + leader + govParty + regionalParl | country, cluster = "parl.tier", data = Model.Data[[i]])
  # Note the different `data` argument here: subset() keeps only national
  # legislators, so this model answers the same question on a narrower sample.
  m4_l[[i]] <- feols(constituency.service ~ TDE + AP + sex + rookie + leader + govParty | country, cluster = "parl.tier", data = subset(Model.Data[[i]], regionalParl == 0))
}

# Pool results across imputations using Rubin's rules
# Use mice::pool to combine results
out1 <- summary(mice::pool(m1_l))
out2 <- summary(mice::pool(m2_l))
out3 <- summary(mice::pool(m3_l))
out4 <- summary(mice::pool(m4_l))

How to read Figure 11.1. One row per survey item, showing how strongly each loads on the underlying constituency-service factor, with uncertainty around each estimate.

This figure validates the measure, not the argument. Large, similarly-signed loadings across the nine items mean they really do reflect a single common trait, which is what licenses collapsing them into one score. An item near zero would be one the factor does not explain — worth noticing, because it says that activity is not part of what the other eight have in common.

Figure 11.1: Constituency Service Factor Loadings

We consider the posterior distribution of factor loadings as a tool to confirm the general appropriateness of our modeling choices. In general, we would expect these not only to be positive, but also clearly bounded away from 0. Outcomes with large factor loadings contribute more heavily to constituency service scores.

# Set bottom margin large enough to display rotated x-axis labels
par(mar = c(7, 4, 0, 0))

# Create empty plotting canvas
plot(c(1, 9), c(0, 1.5),
  xlab = "", ylab = "Factor loadings",
  axes = F, type = "n"
)

# Add y-axis ticks (rotated)
axis(2, las = 1)
# Add x-axis tick positions, but hide default labels
axis(1, at = 1:9, labels = FALSE)

# Add rotated item labels below x-axis
text(seq(1, 9, by = 1), par("usr")[3] - 0.2, labels = colnames(ordFA.loadings), srt = 45, pos = 1, xpd = TRUE)

# Plot posterior means of loadings (dots)
points(xy.coords(1:9, colMeans(ordFA.loadings)), pch = 19)

# Add 95% credible intervals as vertical lines
segments(
  x0 = 1:9, x1 = 1:9,
  y0 = apply(ordFA.loadings, 2, quantile, 0.025),
  y1 = apply(ordFA.loadings, 2, quantile, 0.975)
)

# dev.print(pdf        # copies the plot to a the PDF file
#          , width = 9, height = 7
#          , "factor-loadings-constituency-service.pdf")

The plot displays the factor loadings from the Bayesian ordinal factor model, with credible intervals reflecting uncertainty. All nine items load positively and significantly on the latent factor, confirming that each provides meaningful information about MPs’ engagement in constituency service. As described in the chapter, loadings for items like meetgreet, localmedia, and publicize are close to or above 1, indicating that these behaviors are particularly indicative of high constituency service scores. By contrast, ceremonial or low-visibility behaviors like weddings, meetings, and lectures contribute somewhat less, but still meaningfully. This provides strong support for using a unidimensional latent trait model to summarize MPs’ constituency service activities.

Possible failure point, same family as the chapter 9 error. rubin_out() builds its pooled variance-covariance matrix from vcov() applied to feols models. Recent versions of fixest return that as a classed fixest_vcov object rather than a plain matrix, and simulation functions that inspect the matrix can fail to dispatch on it — the error in chapter 9 read “no applicable method for ‘isSymmetric’ applied to an object of class fixest_vcov”.

Whether it bites here depends on mixtools::rmvnorm(), which is a different implementation from the mvtnorm version that failed in chapter 9. If knitting stops at the Figure 11.3 blocks, the fix is the same one-word change used there — wrap the matrix in unclass():

sim.coefs.4 <- mixtools::rmvnorm(100, correct_m4$betas, unclass(correct_m4$vcov))

It strips a class label and changes no numbers. Left unmodified for now, since it may not be needed.

How to read Table 11.2. Four columns matching the four models above, getting stricter left to right, with the last restricted to national legislators.

Read the AP row across. The chapter’s claim is that higher personal vote-seeking incentives go with more constituency service — but check the sign convention of the factor score before deciding which direction supports it, since the underlying survey items run from most to least frequent.

Then read the TDE row as a placebo: the argument concerns intraparty incentives, so the interparty measure should not be doing the work. Consistency between columns (3) and (4) is the check that the result does not hinge on how regional legislators are handled.

Table 11.2: Association between AP and Constituency Service

To assess the extent to which electoral rules and legislator characteristics explain variation in constituency service, we estimate a series of regression models. Each model uses constituency service factor scores as the dependent variable and includes country fixed effects with clustered standard errors at the parliamentary tier level. AP and TDE are simulated using gradient-boosting machine models based on electoral rules. Results are pooled across five imputations using Rubin’s rules and presented in Table 11.2.

# We need to build a table for the book
model_list <- list(out1, out2, out3, out4)
out_estimates <- sapply(model_list, extract_Estimates, simplify = FALSE)

# Create table
mytable <- merge(out_estimates[[1]], out_estimates[[2]], by = "variable", all = TRUE)
for (i in c(3:4)) mytable <- merge(mytable, out_estimates[[i]], by = "variable", all = TRUE)

# Reorder variables
mytable <- mytable[c(1, 7, 6, 5, 3, 2, 4), ]

# Rename rows
row.names(mytable) <- mytable$variable
mytable$variable <- NULL

# Calculate TSS
TSS <- sum((Model.Data[[1]]$constituency.service - mean(Model.Data[[1]]$constituency.service))^2)
TSS_2 <- sum((na.omit(Model.Data[[1]][, c("constituency.service", "AP", "TDE", "sex", "rookie", "leader", "govParty", "regionalParl")])$constituency.service - mean(na.omit(Model.Data[[1]][, c("constituency.service", "AP", "TDE", "sex", "rookie", "leader", "govParty", "regionalParl")])$constituency.service))^2)
TSS_3 <- sum((na.omit(Model.Data[[1]][Model.Data[[1]]$regionalParl == 0, c("constituency.service", "AP", "TDE", "sex", "rookie", "leader", "govParty", "regionalParl")])$constituency.service - mean(na.omit(Model.Data[[1]][Model.Data[[1]]$regionalParl == 0, c("constituency.service", "AP", "TDE", "sex", "rookie", "leader", "govParty", "regionalParl")])$constituency.service))^2)


# Calculate R2
r2_m1 <- round(Reduce("+", lapply(m1_l, function(x) (1 - (sum(x$residuals^2) / TSS)))) / length(m1_l), 3)
r2_m2 <- round(Reduce("+", lapply(m2_l, function(x) (1 - (sum(x$residuals^2) / TSS)))) / length(m2_l), 3)
r2_m3 <- round(Reduce("+", lapply(m3_l, function(x) (1 - (sum(x$residuals^2) / TSS_2)))) / length(m3_l), 3)
r2_m4 <- round(Reduce("+", lapply(m4_l, function(x) (1 - (sum(x$residuals^2) / TSS_3)))) / length(m4_l), 3)

# Calculate Adj. Pseudo R2
adjr2_m1 <- round(Reduce("+", lapply(m1_l, function(x) (1 - (((sum(x$residuals^2) / TSS) * (x$nobs - 1)) / (x$nobs - x$nparams - 1))))) / length(m1_l), 3)
adjr2_m2 <- round(Reduce("+", lapply(m2_l, function(x) (1 - (((sum(x$residuals^2) / TSS) * (x$nobs - 1)) / (x$nobs - x$nparams - 1))))) / length(m2_l), 3)
adjr2_m3 <- round(Reduce("+", lapply(m3_l, function(x) (1 - (((sum(x$residuals^2) / TSS_2) * (x$nobs - 1)) / (x$nobs - x$nparams - 1))))) / length(m3_l), 3)
adjr2_m4 <- round(Reduce("+", lapply(m4_l, function(x) (1 - (((sum(x$residuals^2) / TSS_3) * (x$nobs - 1)) / (x$nobs - x$nparams - 1))))) / length(m4_l), 3)

# Add R2, adj.R2, and N
mytable <- rbind(mytable, "Fixed Effects by Country" = c("Yes", "Yes", "Yes", "Yes"))
mytable <- rbind(mytable, "Regional Parliaments Included" = c("Yes", "Yes", "Yes", "No"))
mytable <- rbind(mytable, "R2" = c(r2_m1, r2_m2, r2_m3, r2_m4))
mytable <- rbind(mytable, "Adj R2" = c(adjr2_m1, adjr2_m2, adjr2_m3, adjr2_m4))
mytable <- rbind(mytable, "Observations" = c(
  m1_l[[1]]$nobs,
  m2_l[[1]]$nobs,
  m3_l[[1]]$nobs,
  m4_l[[1]]$nobs
))

# Rename rows
row.names(mytable) <- c(
  "AP",
  "TDE",
  "Sex = Male",
  "Rookie",
  "Leader",
  "Gov. Party Member",
  "Regional MP",
  "Fixed Effects by Country",
  "Regional Parliaments Included",
  "R2",
  "Adj. R2",
  "Observations"
)

# Replace NA for '
mytable[is.na(mytable)] <- ""

# Create table
mytable %>%
  kbl(
    caption = ":  Predictors of Constituency Service",
    col.names = c("(1)", "(2)", "(3)", "(4)"),
    format = "latex",
    linesep = "",
  )

Across all models, Average Personalism (AP) emerges as a statistically significant and substantively important predictor of constituency service. Model 1 includes only AP and country fixed effects; Model 2 adds TDE, which is not statistically significant. Model 3 incorporates legislator-level controls, some of which—such as being male, in a leadership position, or in the governing party—are associated with higher scores. Model 4 restricts the sample to national parliaments only, yielding the largest estimated effect of AP (3.595), suggesting that electoral rules may influence constituency service most strongly at the national level.

Figure 11.3: AP and Constituency Service

Figure 11.3(a): effect of AP on constituency service based on national parliaments (model 4), keeping other predictors constant at their median values

Figure 11.3(a) illustrates the marginal effect of Average Personalism (AP) on predicted constituency service scores, based on Model 4. This model includes only national MPs and controls for legislator characteristics such as gender, rookie status, leadership role, and government party membership. The plot uses posterior simulation to estimate expected constituency service levels and a 95% confidence band across a range of observed AP scores.

# Set margin space for axes
par(mar = c(4, 4, 0, 0))
# Define the range of AP scores for the x-axis
newx <- seq(0.35, 0.6, length.out = 100)

# Create a new dataset with AP varying, holding other predictors at their medians
newdata <- data.frame(
  "TDE" = median(c(
    Model.Data[[1]]$TDE, Model.Data[[2]]$TDE,
    Model.Data[[3]]$TDE, Model.Data[[4]]$TDE,
    Model.Data[[5]]$TDE
  )),
  "AP" = newx,
  "sex" = median(Model.Data[[1]]$sex),
  "rookie" = median(Model.Data[[1]]$rookie),
  "leader" = median(Model.Data[[1]]$leader),
  "govParty" = median(Model.Data[[1]]$govParty)
)

# Extract posterior mean and covariance from pooled Rubin output
correct_m4 <- rubin_out(m4_l)


# Simulate 100 draws from the posterior distribution of the coefficients
# Predictions based on robust standard errors
# Each draw is one plausible version of "what the model might be". Feeding all
# 100 through the prediction step and taking percentiles produces the
# confidence bands in Figure 11.3 -- the same simulation approach used in
# Chapter 9. rubin_out() above supplied the pooled coefficients and the pooled
# variance-covariance matrix that this draws from.
# set.seed() fixes the draws so the figure reproduces exactly.
set.seed(123)
sim.coefs.4 <- mixtools::rmvnorm(100, correct_m4$betas, correct_m4$vcov)

# Predict expected values across all AP values
preds.4 <- c()
for (i in 1:100) {
  expected.values.4 <- sim.coefs.4 %*% t(newdata[i, ])
  preds.tmp <- quantile(expected.values.4, prob = c(0.025, 0.5, 0.975))
  preds.4 <- rbind(preds.4, preds.tmp)
}
colnames(preds.4) <- c("lwr", "fit", "upr")


# Plot the marginal effect of AP on constituency service
plot(
  x = c(min(newx), max(newx)), y = c(0, 1),
  xlab = "", ylab = "", type = "n", bty = "n",
  xlim = c(0.35, 0.6), ylim = c(0, 5)
)

# Axis labels
mtext(side = 1, line = 2, text = "Average Personalism")
mtext(side = 2, line = 2, text = "Expected Constituency Service Score")

# Shaded 95% credible interval
polygon(
  x = c(newx, rev(newx)),
  y = c(preds.4[, "lwr"], rev(preds.4[, "upr"])),
  col = adjustcolor("grey", alpha.f = 0.95),
  border = NA
)

# Fitted line
lines(x = newx, y = preds.4[, "fit"], lty = 1, col = "white", lwd = 3)

# dev.print(pdf        # copies the plot to a the PDF file
# , width = 9, height = 7
# , "constituency-service-vs-ap-national.pdf")

This figure shows that as Average Personalism (AP) increases, the expected level of constituency service reported by national MPs rises. The central white line represents the median prediction, and the shaded band represents the 95% credible interval. Moving from the minimum to the maximum observed value of AP results in a sizable increase in predicted constituency service, reinforcing the claim that personalistic electoral environments incentivize greater direct engagement with constituents

How to read Figure 11.3, panel (a). Expected constituency service across the range of AP, with a simulated confidence band, based on the national-parliament model. Judge it by whether the band at one end of the AP range clears the band at the other — that, not the steepness of the line, is what separates the relationship from noise.

Figure 11.3b: probability of engaging in various constituency service activities at least once every 2 weeks based on minimum and maximum AP scores.

Figure 11.3(b) provides a more granular interpretation of the association between Average Personalism (AP) and constituency service. Rather than summarizing engagement as a single latent score (as in previous figures), this plot decomposes the effect across the nine observed activities from the PARTIREP survey. Using the results from the ordered factor analysis and the predictive model (Model 4), the plot shows the predicted probability that a prototypical MP (male, experienced, non-leader, from a government party) will report engaging in each activity at least once every two weeks, under two scenarios: low AP (min value in the data) and high AP (max value).

# Extract factor loadings, thresholds, and gamma values from the ordinal factor analysis
# These define how each observed item relates to the latent trait

ordFA.loadings <- ordFactAnal[, grep("Lambda", colnames(ordFactAnal))]
ordFA.lambda.1 <- ordFA.loadings[, grep("\\.1", colnames(ordFA.loadings))]
ordFA.lambda.2 <- ordFA.loadings[, grep("\\.2", colnames(ordFA.loadings))]
ordFA.gamma <- ordFactAnal[, grep("gamma", colnames(ordFactAnal))]
ordFA.gamma.var1 <- ordFA.gamma[, grep("\\_1", colnames(ordFA.gamma))]
ordFA.gamma.var2 <- ordFA.gamma[, grep("\\_2", colnames(ordFA.gamma))]
ordFA.gamma.var3 <- ordFA.gamma[, grep("\\_3", colnames(ordFA.gamma))]
ordFA.gamma.var4 <- ordFA.gamma[, grep("\\_4", colnames(ordFA.gamma))]
ordFA.gamma.var5 <- ordFA.gamma[, grep("\\_5", colnames(ordFA.gamma))]
ordFA.gamma.var6 <- ordFA.gamma[, grep("\\_6", colnames(ordFA.gamma))]
ordFA.gamma.var7 <- ordFA.gamma[, grep("\\_7", colnames(ordFA.gamma))]
ordFA.gamma.var8 <- ordFA.gamma[, grep("\\_8", colnames(ordFA.gamma))]
ordFA.gamma.var9 <- ordFA.gamma[, grep("\\_9", colnames(ordFA.gamma))]

ordFA.gamma.list <- list(
  ordFA.gamma.var1,
  ordFA.gamma.var2,
  ordFA.gamma.var3,
  ordFA.gamma.var4,
  ordFA.gamma.var5,
  ordFA.gamma.var6,
  ordFA.gamma.var7,
  ordFA.gamma.var8,
  ordFA.gamma.var9
)

# We choose min and max sample values of AP (pers.hat) to illustrate effects
ap <- c(
  min(c(
    Model.Data[[1]]$AP, Model.Data[[2]]$AP,
    Model.Data[[3]]$AP, Model.Data[[4]]$AP,
    Model.Data[[5]]$AP
  )),
  max(c(
    Model.Data[[1]]$AP, Model.Data[[2]]$AP,
    Model.Data[[3]]$AP, Model.Data[[4]]$AP,
    Model.Data[[5]]$AP
  ))
)


# Set constant values of intercept, tde, ap, male, rookie, leader
# govParty
set.values <- cbind(
  c(0.631, ap[1], 1, 0, 0, 1),
  c(0.631, ap[2], 1, 0, 0, 1)
)


# Produce expected constituency scores
# (one column for each of the interesting values of AP)
expected.constituency.score <- sim.coefs.4 %*% set.values

# Set relevant parameters for outcome variables (choose gamma)
gamma.var.num <- c("\\_1", "\\_2", "\\_3", "\\_5", "\\_8", "\\_4", "\\_6", "\\_7", "\\_9") # Same order of presentation of variables

probs.lo.ap <- probs.hi.ap <- matrix(NA, ncol = 3, nrow = 9)
for (j in 1:9) { # We only need the third column of Gamma, because
  # we'll look at how the probability of choosing                        # 5 or 6 changes with changes in AP
  tmp.lo <- tmp.hi <- c()
  for (i in 1:100) {
    tmp.1 <- 1 - pnorm(ordFA.gamma.list[[j]][, 3] - ordFA.lambda.1[, grep(gamma.var.num[j], colnames(ordFA.lambda.1))] - ordFA.lambda.2[, grep(gamma.var.num[j], colnames(ordFA.lambda.2))] * expected.constituency.score[i, 1])
    tmp.2 <- 1 - pnorm(ordFA.gamma.list[[j]][, 3] - ordFA.lambda.1[, grep(gamma.var.num[j], colnames(ordFA.lambda.1))] - ordFA.lambda.2[, grep(gamma.var.num[j], colnames(ordFA.lambda.2))] * expected.constituency.score[i, 2])
    tmp.lo <- c(tmp.lo, tmp.1)
    tmp.hi <- c(tmp.hi, tmp.2)
  }
  probs.lo.ap[j, ] <- quantile(tmp.lo, probs = c(0.25, 0.5, 0.75)) # 50% CI
  probs.hi.ap[j, ] <- quantile(tmp.hi, probs = c(0.25, 0.5, 0.75))
}

# Increase in probability of engaging at least once every two weeks in each activity
mean(probs.hi.ap[, 2])

## [1] 0.684

mean(probs.lo.ap[, 2])

## [1] 0.529

# Plot of probabilities of engaging at least once every two weeks
par(mar = c(5, 4, 0, 1))
plot(c(1, 9), c(0, 1),
  xlab = "", ylab = "Probability of engaging at least once every two weeks",
  axes = F, type = "n"
)

# Add axes
axis(2, las = 1)
axis(1, at = 1:9, labels = FALSE)
text(seq(1, 9, by = 1), par("usr")[3] - 0.08,
  labels = c(
    "weddings", "meetings", "lectures",
    "surgery", "meetgreet", "newslttr",
    "advertise", "publicize", "localmedia"
  ),
  srt = 45, pos = 1, xpd = TRUE
)

# Plot median probabilities and 50% intervals (black = high AP, gray = low AP)
points(xy.coords(c(1:9) - 0.1, probs.lo.ap[, 2]), pch = 19, col = "gray")
points(xy.coords(c(1:9) + 0.1, probs.hi.ap[, 2]), pch = 19, col = "black")
segments(
  x0 = c(1:9) - 0.1, x1 = c(1:9) - 0.1,
  y0 = probs.lo.ap[, 1],
  y1 = probs.lo.ap[, 3],
  col = "gray"
)
segments(
  x0 = c(1:9) + 0.1, x1 = c(1:9) + 0.1,
  y0 = probs.hi.ap[, 1],
  y1 = probs.hi.ap[, 3],
  col = "black"
)
# Add legend
legend("topleft",
  bty = "n", pch = 19, col = c("gray", "black"),
  legend = c("Low AP", "High AP")
)

# dev.print(pdf        # copies the plot to a the PDF file
#           , width = 9, height = 7
#           , "predicted-constituency-service-activities.pdf")

This figure shows that across all nine indicators of constituency service, the probability of frequent engagement increases with AP. Under low AP electoral environments, the average probability of engaging in any activity at least once every two weeks is around 0.52; under high AP, this rises to 0.68. The largest increases in probability are observed in activities like: meetgree (meeting local businesses and groups), localmedia (media presence), publicize (credit-claiming efforts). These suggest that the personal vote is strongly tied to visibility and outreach efforts. More ceremonial or reactive forms of service, such as weddings and lectures, see more modest increases.

Average correlation

# AP and TDE
mean(sapply(1:5, function(i) cor(Model.Data[[i]]$AP, Model.Data[[i]]$TDE)))

## [1] 0.0928

mean(sapply(1:5, function(i) cor(Model.Data[[i]]$constituency.service, Model.Data[[i]]$TDE)))

## [1] 0.125

These low positive correlations confirm the discussion in the chapter: there is some overlap between Average Personalism (AP) and Total Duvergerian Effect (TDE), likely because SMDP and CLPR systems co-occur frequently in the data. However, the relationship is weak, and TDE is not a robust predictor of constituency service. This supports its inclusion as a control, not a theoretically motivated covariate.

How to read Figure 11.3, panel (b). Panel (a) worked in units of the factor score, which has no natural meaning. This panel translates the same result back into the individual activities the score was built from: the probability that a legislator engages in each one, at low versus high AP.

That translation is what makes the finding legible. “A 0.3 change in a latent factor score” means little; “legislators under high-AP rules are more likely to hold surgeries and attend local events” is a claim a reader can weigh.

Ambition as a potential driver of results

A rival explanation, tested. Perhaps legislators who do more constituency service are simply the more ambitious ones, and ambition happens to be distributed differently across electoral systems — in which case the electoral rules would not be doing the causal work. This short section checks that possibility rather than leaving it to the reader.

This block compares AP distribution and political ambition levels between national and regional MPs.

# Distribution of AP
# Among national MPs
summary(Model.Data[[1]]$AP[Model.Data[[1]]$regionalParl == 0])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.354   0.391   0.446   0.438   0.480   0.575

sd(Model.Data[[1]]$AP[Model.Data[[1]]$regionalParl == 0], na.rm = T)

## [1] 0.0567

# Among regional MPs
summary(Model.Data[[1]]$AP[Model.Data[[1]]$regionalParl == 1])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.358   0.392   0.470   0.451   0.482   0.537

sd(Model.Data[[1]]$AP[Model.Data[[1]]$regionalParl == 1], na.rm = T)

## [1] 0.0513

# national MPs
table(Model.Data[[1]]$ambition[Model.Data[[1]]$regionalParl == 0])

## 
##   0   1 
## 318 529

mean(Model.Data[[1]]$ambition[Model.Data[[1]]$regionalParl == 0], na.rm = T)

## [1] 0.625

# regional MPs
table(Model.Data[[1]]$ambition[Model.Data[[1]]$regionalParl == 1])

## 
##   0   1 
## 486 668

mean(Model.Data[[1]]$ambition[Model.Data[[1]]$regionalParl == 1], na.rm = T)

## [1] 0.579

AP is slightly higher among regional MPs than national MPs (mean of 0.451 vs. 0.438), contrary to the idea that national MPs operate in more personalistic systems. Standard deviations and ranges are similar, indicating comparable variation across groups. Political ambition is marginally higher among national MPs (0.625 vs. 0.579), consistent with the idea that national legislators might be more career-oriented—but the difference is small.

Figure 11.4: Constituency Service, District Magnitude and Ballot Type

Figure 11.4 explores how the effect of electoral personalism on constituency service varies across combinations of ballot type (Closed vs. Open List PR) and district magnitude (M). The figure uses the simulation-based AP estimates from the GBM models to compute predicted constituency service scores across six values of M (2–100), holding other institutional and individual-level predictors constant. This approach tests Carey and Shugart’s (1995) conjecture that incentives to cultivate a personal vote increase with M under open lists, and are generally lower under closed lists.

# Define electoral rule scenarios: Hare quota, one vote, party-level pooling
# Vary ballot type and M
xAPclosed <- data.frame(
  formula = "hare",
  ballot_type = c("closed"),
  new.nvotes = "One",
  pool_level = "party",
  threshold = 0,
  M = c(2, 10, 15, 25, 50, 100)
)

# Predict AP scores using GBM models for each ballot type and M combination
xAPopen <- data.frame(
  formula = "hare",
  ballot_type = c("open"),
  new.nvotes = "One",
  pool_level = "party",
  threshold = 0,
  M = c(2, 10, 15, 25, 50, 100)
)
xAPclosed$formula <- as.factor(xAPclosed$formula)
xAPclosed$ballot_type <- as.factor(xAPclosed$ballot_type)
xAPclosed$new.nvotes <- as.factor(xAPclosed$new.nvotes)
xAPclosed$pool_level <- as.factor(xAPclosed$pool_level)
xAPopen$formula <- as.factor(xAPopen$formula)
xAPopen$ballot_type <- as.factor(xAPopen$ballot_type)
xAPopen$new.nvotes <- as.factor(xAPopen$new.nvotes)
xAPopen$pool_level <- as.factor(xAPopen$pool_level)

# Combine both
xAP <- rbind(xAPclosed, xAPopen)

# Get TDE, these create lists
ap_hat1 <- predict(optimalGBMIntra[[1]], xAP, n.trees = optimalGBMIntra[[1]]$n.trees)
ap_hat2 <- predict(optimalGBMIntra[[2]], xAP, n.trees = optimalGBMIntra[[2]]$n.trees)
ap_hat3 <- predict(optimalGBMIntra[[3]], xAP, n.trees = optimalGBMIntra[[3]]$n.trees)
ap_hat4 <- predict(optimalGBMIntra[[4]], xAP, n.trees = optimalGBMIntra[[4]]$n.trees)
ap_hat5 <- predict(optimalGBMIntra[[5]], xAP, n.trees = optimalGBMIntra[[5]]$n.trees)

# Add TDE and AP
xAP$ap <- colMeans(rbind(ap_hat1, ap_hat2, ap_hat3, ap_hat4, ap_hat5))


# Create two data frames for predicting constituency service under each scenario
outClose.data <- data.frame(
  "TDE" = median(c(
    Model.Data[[1]]$TDE, Model.Data[[2]]$TDE,
    Model.Data[[3]]$TDE, Model.Data[[4]]$TDE,
    Model.Data[[5]]$TDE
  )),
  "AP" = xAP$ap[xAP$ballot_type == "closed"],
  "sex" = median(Model.Data[[1]]$sex),
  "rookie" = median(Model.Data[[1]]$rookie),
  "leader" = median(Model.Data[[1]]$leader),
  "govParty" = median(Model.Data[[1]]$govParty)
)

outOpen.data <- data.frame(
  "TDE" = median(c(
    Model.Data[[1]]$TDE, Model.Data[[2]]$TDE,
    Model.Data[[3]]$TDE, Model.Data[[4]]$TDE,
    Model.Data[[5]]$TDE
  )),
  "AP" = xAP$ap[xAP$ballot_type == "open"],
  "sex" = median(Model.Data[[1]]$sex),
  "rookie" = median(Model.Data[[1]]$rookie),
  "leader" = median(Model.Data[[1]]$leader),
  "govParty" = median(Model.Data[[1]]$govParty)
)

outOpen <- c()
outClose <- c()

# Simulate predicted constituency scores using posterior coefficients
for (i in 1:6) {
  exp.open <- sim.coefs.4 %*% t(outOpen.data[i, ])
  exp.close <- sim.coefs.4 %*% t(outClose.data[i, ])
  preds.open <- quantile(exp.open, prob = c(0.25, 0.5, 0.75))
  preds.close <- quantile(exp.close, prob = c(0.25, 0.5, 0.75))
  outOpen <- rbind(outOpen, preds.open)
  outClose <- rbind(outClose, preds.close)
}
colnames(outClose) <- colnames(outOpen) <- c("lwr", "fit", "upr")

M <- c(2, 10, 15, 25, 50, 100)

# Plot: Expected constituency service score vs. District Magnitude
plot(
  x = c(1, 6), y = c(1, 3),
  xlab = "", ylab = "", type = "n",
  bty = "n", axes = F
)
axis(2)
axis(1, at = c(1:6), labels = M)
mtext(side = 1, line = 2.5, text = "District magnitude", cex = 1.5)
mtext(side = 2, line = 2.5, text = "Expected constituency service score", cex = 1.5)
segments(
  x0 = c(1:6) - 0.1, x1 = c(1:6) - 0.1,
  y0 = outClose[, 1], y1 = outClose[, 3],
  lwd = 3, col = "gray"
)
segments(
  x0 = c(1:6) + 0.1, x1 = c(1:6) + 0.1,
  y0 = outOpen[, 1], y1 = outOpen[, 3],
  lwd = 3, col = "black"
)
points(xy.coords(c(1:6) - 0.1, outClose[, 2]), pch = 19, col = "gray")
points(xy.coords(c(1:6) + 0.1, outOpen[, 2]), pch = 19, col = "black")
legend("topright",
  pch = c(19, 19), col = c("gray", "black"),
  c("CLPR", "OLPR"),
  bty = "n"
)

# dev.print(pdf        # copies the plot to a the PDF file
# , width = 9, height = 7
# , "constituency-service-closed-vs-open-rule.pdf")

This figure shows that, across all levels of district magnitude, Open List PR (OLPR) systems consistently generate higher predicted constituency service scores than Closed List PR (CLPR) systems. This supports the theoretical expectation that ballot structure—particularly the degree to which voters can influence which candidates are elected—affects the incentives legislators have to cultivate a personal vote. However, the results do not support Carey and Shugart’s (1995) expectation that constituency service should increase as M increases under OLPR. Instead: Predicted scores in OLPR systems appear fairly stable across district sizes, with only a slight decline after M = 25.

Under CLPR, predicted scores are consistently lower and also flat or declining with M. This finding suggests that while ballot openness matters, district magnitude does not systematically increase personal vote incentives beyond a moderate threshold. This contradicts classic expectations, reinforcing the book’s broader claim that personal vote incentives are complex and not strictly additive across institutional dimensions

How to read Figure 11.4. District magnitude on the horizontal axis, one line per ballot type. The vertical gap between the lines is the effect of ballot structure; whether that gap widens or narrows across magnitude is the interaction.

This is the chapter’s engagement with the classic Carey and Shugart expectation that magnitude amplifies the personalising effect of open ballots. Roughly parallel lines would say ballot type matters but magnitude does not condition it.

Carey and Shugart Model

The literature’s own specification. The models below use observed district magnitude and ballot type directly, interacted, rather than the GBM-derived AP score used everywhere else. That is deliberate: it tests Carey and Shugart’s argument on its own terms, so the comparison with the AP-based results is a like-for-like one about which measure better explains constituency service.

To further assess whether electoral systems generate stronger incentives for constituency service as district magnitude increases—especially under open lists—we fit two linear interaction models. These regress constituency service scores on the interaction between district magnitude (M) and ballot type, with and without additional controls. These models provide a direct statistical test of Carey and Shugart’s (1995) prediction that larger districts under open-list systems (OLPR) will generate more personal vote-seeking behavior. Both models include country fixed effects, and robust standard errors are used to account for heteroskedasticity.

pdf("~/Downloads/ConstituencyService_CareyShugart.pdf")

# Based on Model 1
# Model 1: Only the interaction between district magnitude and ballot type

c_and_s_m1 <- feols(constituency.service ~ M * ballot_type | country, se = "hetero", data = subset(Model.Data[[1]]))

# Plot interaction: baseline model
sjPlot::plot_model(c_and_s_m1,
  type = "int", terms = c("M", "ballot_type"),
  axis.title = c("M", "Constituency Service"),
  title = "Expected Value of Constituency Service - With Country FE, No Controls",
  legend.title = "Ballot Type"
)
summary(c_and_s_m1)

## OLS estimation, Dep. Var.: constituency.service
## Observations: 2,286
## Fixed-effects: country: 15
## Standard-errors: Heteroskedasticity-robust 
##                        Estimate Std. Error t value   Pr(>|t|)    
## M                      0.001099    0.00114  0.9664 0.33393613    
## ballot_typeflexible    0.334927    0.09745  3.4368 0.00059931 ***
## ballot_typeopen        0.130402    0.13523  0.9643 0.33498844    
## M:ballot_typeflexible -0.009640    0.00234 -4.1132 0.00004041 ***
## M:ballot_typeopen     -0.000183    0.00217 -0.0843 0.93281399    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.703318     Adj. R2: 0.382673
##                  Within R2: 0.011277

# Based on Model 3
# Model 3: Full specification with controls
c_and_s_m3 <- feols(constituency.service ~ M * ballot_type + sex + rookie + leader + govParty + regionalParl | country, se = "hetero", data = subset(Model.Data[[1]]))

# Plot interaction: full model
sjPlot::plot_model(c_and_s_m3,
  type = "int", terms = c("M", "ballot_type"),
  axis.title = c("M", "Constituency Service"),
  title = "Expected Value of Constituency Service -  With Country FE and Controls",
  legend.title = "Ballot Type"
)
summary(c_and_s_m3)

## OLS estimation, Dep. Var.: constituency.service
## Observations: 2,240
## Fixed-effects: country: 15
## Standard-errors: Heteroskedasticity-robust 
##                       Estimate Std. Error t value   Pr(>|t|)    
## M                      0.00205    0.00135   1.527 1.2690e-01    
## ballot_typeflexible    0.30455    0.09598   3.173 1.5291e-03 ** 
## ballot_typeopen        0.15904    0.14028   1.134 2.5705e-01    
## sex                    0.19015    0.03233   5.881 4.6843e-09 ***
## rookie                 0.03599    0.03061   1.176 2.3988e-01    
## leader                 0.10042    0.04441   2.261 2.3845e-02 *  
## govParty               0.10335    0.03178   3.252 1.1638e-03 ** 
## regionalParl          -0.06317    0.04763  -1.326 1.8495e-01    
## M:ballot_typeflexible -0.01014    0.00235  -4.318 1.6460e-05 ***
## M:ballot_typeopen     -0.00127    0.00230  -0.551 5.8170e-01    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.696522     Adj. R2: 0.383345
##                  Within R2: 0.034789

dev.off()

## quartz_off_screen 
##                 2

The results from both interaction models offer no support for Carey and Shugart’s (1995) hypothesis that larger district magnitudes increase personal vote incentives under open-list systems. In both the baseline model with country fixed effects and the fully controlled model, the interaction between district magnitude (M) and open-list ballot structures is statistically insignificant, while the interaction with flexible lists is significantly negative. This suggests that in flexible-list systems, constituency service actually declines as M increases. Moreover, the effect of M alone is not significant in either model, and the primary determinant of constituency service remains the ballot type itself. These findings align with the simulation-based results from Figure 11.4 and reinforce the broader conclusion that ballot openness—not district size—is the dominant institutional driver of personal vote-seeking behavior.

Moving Forward

Chapter 11 shows that electoral system incentives, measured by Average Personalism (AP), strongly influence legislators’ engagement in constituency service. Using survey data from the PARTIREP project, the chapter finds that MPs in candidate-centered systems report more frequent constituency service, highlighting how intraparty competition shapes politicians’ efforts to build personal support. This relationship holds after controlling for individual, party, and institutional factors and is stronger among national legislators than regional ones. The findings underscore the importance of electoral rules not just for who gets elected, but for how representatives behave once in office.

The next chapter shifts focus to legislative institutions, exploring how intraparty incentives influence the organization and functioning of legislatures. Together, these chapters deepen our understanding of how electoral design shapes political behavior both during campaigns and in governance.

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
coda	0.19-4.1
dplyr	1.1.4
fixest	0.14.1
forcats	1.0.0
gbm	2.2.2
ggplot2	3.5.2
kableExtra	1.4.1
labelled	2.16.0
lubridate	1.9.4
MASS	7.3-65
MCMCpack	1.7-1
mice	3.19.0
mixtools	2.0.0.1
psych	2.6.5
purrr	1.2.2
readr	2.1.5
readxl	1.4.5
sandwich	3.1-1
sjPlot	2.9.0
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.

11.Constituency Service