Waking Up the Golden Dawn
Does Exposure to the Refugees Increase Support for Extreme-Right Parties?
The 2015 refugee crisis was the largest displacement in Europe since World War II, yet its political consequences remain contested. This article asks whether sudden exposure to refugees increased support for Golden Dawn, Greece’s extreme-right party. Leveraging a natural experiment created by geography, I analyze four elections (2012–2015) with difference-in-differences and event-study models. Results show that even brief exposure—without sustained contact or competition—raised Golden Dawn’s vote share by about two percentage points (44%). The findings highlight the “flash potential” of refugee crises to reshape politics and challenge established theories of contact and conflict.
This document is designed as a didactic model for how to structure and write a quantitative research paper using Quarto, LaTeX, Jabref, and credible identification strategies such as difference-in-differences and event-studies. It demonstrates best practices in framing a research question, engaging with the literature, building theory, and presenting empirical evidence. For reasons of brevity, the paper reproduces only a subset of the original results, serving primarily as a teaching and training resource. The full study and findings are available in the published article:
Dinas E, Matakos K, Xefteris D, Hangartner D. Waking Up the Golden Dawn: Does Exposure to the Refugee Crisis Increase Support for Extreme-Right Parties? Political Analysis. 2019;27(2):244-254. doi:10.1017/pan.2018.48
1 Introduction
Context
The refugee crisis is the largest displacement of people in Europe since the aftermath of World War II. Since 2015, millions of asylum seekers—primarily Syrians—have crossed the Mediterranean into Europe (UNHCR 2017). Host countries have faced humanitarian, economic, and political challenges (Foged and Peri 2016), including protests and violence against refugees (Editorial Board 2015). These developments have intensified debates over how exposure to refugees shapes electoral support for extreme-right parties.
Research Question
This article investigates whether exposure to the refugee crisis increased support for the extreme-right party Golden Dawn in Greece. Specifically: Did municipalities in Greece that received sudden refugee arrivals in 2015 experience higher vote shares for Golden Dawn than comparable municipalities that did not?
Importance
Answering this question matters both empirically and theoretically. While a large literature shows that labor migration often fuels anti-immigrant parties (Barone et al. 2016; Mendez and Cutillas 2014; Halla, Wagner, and Zweimüller 2015; Brunner and Kuhn 2014; Becker and Fetzer 2016), less is known about refugees. Refugees differ in motives and in the humanitarian concerns they evoke (Bansak, Hainmueller, and Hangartner 2016). Recent studies in Austria and Denmark reach mixed conclusions (Steinmayr 2016; Dustmann, Vasiljeva, and Damm 2016), leaving a gap in our understanding—especially for frontline states like Greece, which received more than half of all arrivals in 2015 (UNHCR 2016).
Methods and Data
The analysis leverages a natural experiment created by geography: some Aegean islands near the Turkish coast received massive inflows of refugees, while nearby islands with similar characteristics did not. This quasi-random exposure enables causal identification. The geographic distribution of arrivals is shown in Figure 1, which highlights the stark contrast between municipalities that were exposed and those that were not. Electoral outcomes across four Greek parliamentary elections (2012–2015) are analyzed using difference-in-differences (Bertrand, Duflo, and Mullainathan 2004).
Show the code
# ------------------------------------------------------------
# Goal: Read election + map data for Greece, aggregate to the
# municipality-year level, compute vote shares/turnout,
# merge with geographies, and plot refugee arrivals per capita.
# This version adds step-by-step comments for R beginners.
# ------------------------------------------------------------
# --- 0) Housekeeping ---------------------------------------------------------
# Remove all objects from memory so the script starts fresh.
rm(list = ls())
# --- 1) Load packages --------------------------------------------------------
# "haven" reads .dta (Stata) files. "dplyr" is for data wrangling.
# "sf" handles spatial (map) data. "ggplot2" makes plots.
# "here" helps build file paths that work on any computer.
# "ggpubr" arranges multiple ggplots in a grid.
library(haven)
library(dplyr)
library(sf)
sf::sf_use_s2(FALSE)  # turn off spherical geometry to avoid certain topology errors
library(ggplot2)
library(here)
library(ggpubr)
# --- 2) Read data from the /data_final folder -------------------------------
# NOTE: `here("data_final", "file.ext")` builds a path like
# "<your-project-root>/data_final/file.ext". Make sure your working
# directory is the project root (where your .Rproj lives).
countries          <- read_sf(here("data_final", "countries.shp"))
gr_muni            <- read_sf(here("data_final", "gr_muni.shp"))
gr_lemnos_island   <- read_sf(here("data_final", "gr_lemnos_island.shp"))
gr_limnos_island   <- read_sf(here("data_final", "gr_limnos_island.shp"))
all_data           <- read_dta(here("data_final", "dinas_data.dta"), encoding = "UTF-8")
# TIP: If you get errors here, check that files exist and that all
# shapefiles share the same CRS (coordinate reference system).
# Use st_crs(object) to inspect and st_transform(...) to convert if needed.
# --- 3) Sort + group + compute totals at municipality-year level ------------
# We first sort the raw data (not strictly necessary, but helpful to read),
# then group by municipality & year so sums are within each municipality-year.
# The column names below (e.g., `valid`, `gd`, `nd`, etc.) must exist in your data.
all_data <- all_data %>%
  arrange(township, municipality, perfecture, year) %>%   # keep original order logic
  group_by(municipality, year) %>%                        # group for within-muni-year sums
  mutate(
    sumall    = sum(valid,      na.rm = TRUE),  # total valid votes
    sumgd     = sum(gd,         na.rm = TRUE),  # Golden Dawn votes (example)
    sumnd     = sum(nd,         na.rm = TRUE),  # New Democracy votes (example)
    sumsyriza = sum(syriza,     na.rm = TRUE),
    sumpsk    = sum(pasok,      na.rm = TRUE),  # PASOK
    sumkke    = sum(kke,        na.rm = TRUE),
    sumanel   = sum(anel,       na.rm = TRUE),
    sumreg    = sum(registered, na.rm = TRUE)   # number registered to vote
  )
# --- 4) Create vote shares + turnout ----------------------------------------
# Shares = party_total / all valid votes. Turnout = valid / registered.
# Watch for division by zero; NA/NaN can occur if `sumall` or `sumreg` is 0.
all_data <- all_data %>%
  mutate(
    gdvote      = sumgd     / sumall,
    ndvote      = sumnd     / sumall,
    syrizavote  = sumsyriza / sumall,
    pasokvote   = sumpsk    / sumall,
    kkevote     = sumkke    / sumall,
    anelvote    = sumanel   / sumall,
    turnout     = sumall    / sumreg
  )
# --- 5) Keep a single row per municipality-year -----------------------------
# The raw file may have multiple rows per municipality-year (e.g., by township).
# We create a within-group row number and keep the first one (any will do
# once we've already computed muni-year sums above).
datamuni <- all_data %>%
  group_by(geo_muni_id, year) %>%  # use the stable geocode id for join later
  mutate(id = dplyr::row_number()) %>%
  filter(id == 1) %>%              # keep one row per geo_muni_id-year
  arrange(geo_muni_id, year) %>%
  group_by(geo_muni_id) %>%
  mutate(id2 = dplyr::row_number())  # a running index within each municipality
# --- 6) Focus on the 2016 election ------------------------------------------
datamuni_2016 <- subset(datamuni, year == 2016)
# --- 7) Merge tabular data to the municipality shapes ------------------------
# left_join keeps all geometries from gr_muni and adds the 2016 data columns.
# Keys: `ge_mn_d` in the shapefile == `geo_muni_id` in the data.
greece_merge <- left_join(gr_muni, datamuni_2016, by = c("ge_mn_d" = "geo_muni_id"))
# --- 8) Handle missing + extreme values in arrivals per capita --------------
# Remove rows where `trarrprop` (arrivals per capita) is missing so the map
# has valid values to color. Then cap extreme values to improve readability.
greece_merge <- subset(greece_merge, !is.na(trarrprop))
# Make a working copy for plotting + truncation
greece_merge2 <- greece_merge
# One municipality has a very large value (> 100). Following the paper:
# "We truncate the variable at 5 arrivals per capita (i.e., 5 refugees per
# resident), which corresponds to the 99.7th percentile of the non‑zero distribution."
# Doing this avoids a legend dominated by outliers.
greece_merge2$trarrprop[greece_merge2$trarrprop > 100] <- 5
# --- 9) Build two maps (zoomed-in Aegean and full Greece) -------------------
# NOTE: `geom_sf()` draws map layers. The `fill` aesthetic colors municipalities
# by `trarrprop`. The color scale runs from green (low) to red (high).
# (a) Zoomed map around 24–29 E longitude, 35.8–40.2 N latitude
fig1a1 <- ggplot() +
  geom_sf(data = countries, fill = "grey80") +
  geom_sf(data = greece_merge2, aes(fill = trarrprop)) +
  scale_fill_gradient(
    low = "green", high = "red", na.value = "white",
    name = "Refugee Arrivals\n(p.c.)\nJan–Sep 2015"
  ) +
  # Draw and label Lemnos/Limnos to match the original figure
  geom_sf(data = gr_lemnos_island,  color = "black", fill = NA, linewidth = 0.3) +
  geom_sf(data = gr_limnos_island,  color = "black", fill = NA, linewidth = 0.3) +
  geom_sf_label(data = gr_limnos_island, aes(label = name),
                size = 3, nudge_x = 0.35, nudge_y = 0.35, alpha = 0.2) +
  geom_sf_label(data = gr_lemnos_island, aes(label = name),
                size = 3, nudge_x = -0.35, nudge_y = -0.35, alpha = 0.2) +
  coord_sf(xlim = c(24, 29), ylim = c(35.8, 40.2)) +
  theme_bw() +
  theme(
    legend.title        = element_text(size = 12, face = "bold"),
    legend.text         = element_text(size = 10),
    legend.key          = element_rect(fill = "yellow", color = "grey", size = 1),
    legend.key.size     = unit(1.5, "lines"),
    legend.background   = element_rect(fill = "white", color = "grey"),
    legend.box.background = element_rect(color = "grey"),
    # Position legend inside the plotting area at the top-right corner
    legend.position     = c(1, 1),
    legend.justification = c("right", "top"),
    axis.title.x        = element_blank(),
    axis.title.y        = element_blank()
  )
# (b) Wider map of Greece
fig1a2 <- ggplot() +
  geom_sf(data = countries, fill = "grey80") +
  geom_sf(data = greece_merge2, aes(fill = trarrprop)) +
  scale_fill_gradient(
    low = "green", high = "red", na.value = "white",
    name = "Refugee Arrivals\n(p.c.)\nJan–Sep 2015"
  ) +
  geom_sf(data = gr_lemnos_island, color = "black", fill = NA, linewidth = 0.3) +
  geom_sf(data = gr_limnos_island, color = "black", fill = NA, linewidth = 0.3) +
  # Labels are optional here; commented out to reduce clutter
  # geom_sf_label(data = gr_limnos_island, aes(label = name),  size = 3, nudge_x = 0.35, nudge_y = 0.35, alpha = 0.2) +
  # geom_sf_label(data = gr_lemnos_island, aes(label = name),  size = 3, nudge_x = -0.35, nudge_y = -0.35, alpha = 0.2) +
  coord_sf(xlim = c(20, 29), ylim = c(34.5, 42)) +
  theme_bw() +
  theme(
    legend.title        = element_text(size = 12, face = "bold"),
    legend.text         = element_text(size = 10),
    legend.key          = element_rect(fill = "yellow", color = "grey", size = 1),
    legend.key.size     = unit(1.5, "lines"),
    legend.background   = element_rect(fill = "white", color = "grey"),
    legend.box.background = element_rect(color = "grey"),
    legend.position     = c(1, 1),
    legend.justification = c("right", "top"),
    axis.title.x        = element_blank(),
    axis.title.y        = element_blank()
  )
# --- 10) Arrange the two plots side-by-side ----------------------------------
# `ggarrange()` puts multiple ggplots in a single figure.
# `common.legend = TRUE` puts one shared legend at the bottom.
row1 <- ggarrange(fig1a1, fig1a2, ncol = 2, common.legend = TRUE, legend = "bottom")
# --- 11) Optional: save to file ---------------------------------------------
# ggsave(here("figures", "refugee_arrivals_maps.png"), row1, width = 10, height = 6, dpi = 300)
# Print to the RStudio plotting pane (or your device if using ggsave later)
print(row1)Findings and Contribution
Results show that exposure alone, even without sustained contact or competition for resources, increased support for Golden Dawn by around 2 percentage points—a 44 percent rise relative to its baseline. Refugees’ presence on islands was typically very brief (Capon 2015), making these findings difficult to reconcile with contact theory (Allport 1979) or realistic group conflict theory (Campbell 1965). Instead, they support arguments about the “flash potential” of refugee crises to reshape politics (Sniderman, Hagendoorn, and Prior 2004; Hopkins 2010).
Structure of the Paper
The following sections situate the case in the broader European context, describe the data and identification strategies in detail, present the empirical results, and conclude by discussing their theoretical and political implications.
2 Literature Review
Introduction to the Debate
The relationship between migration and support for extreme-right parties has been widely studied, but the findings differ depending on whether the focus is on labor migration or refugee migration. While much of the literature documents a positive association between immigration and support for far-right parties, fewer studies isolate the specific effects of refugee inflows. This distinction matters because refugees are perceived differently from labor migrants: their arrival is often framed in humanitarian rather than economic terms (Bansak, Hainmueller, and Hangartner 2016).
Labor Migration and Electoral Backlash
A large body of research links labor migration to electoral gains for anti-immigrant parties. Studies across Spain, Austria, Germany, and Switzerland find that increases in immigration correlate with stronger support for the far-right (Barone et al. 2016; Mendez and Cutillas 2014; Halla, Wagner, and Zweimüller 2015; Brunner and Kuhn 2014; Becker and Fetzer 2016). These works often emphasize competition over jobs, housing, and welfare resources, consistent with realistic group conflict theory (Campbell 1965). The underlying claim is that immigration intensifies competition, leading natives to support exclusionary parties.
Refugees as a Distinct Category
Refugee inflows differ from labor migration in both scale and framing. Humanitarian concerns can make host populations more sympathetic (Bansak, Hainmueller, and Hangartner 2016), potentially moderating backlash. Recent work has directly examined refugees’ political effects, though findings remain mixed. For Austria, exposure to refugees reduced support for the far-right, consistent with contact theory, which argues that sustained interaction with out-groups reduces prejudice (Steinmayr 2016; Allport 1979). By contrast, in Denmark, exposure to refugees increased hostility and electoral support for the far-right (Dustmann, Vasiljeva, and Damm 2016). These contradictory results suggest that context and mechanisms of exposure matter: whether refugees are seen as temporary transients or potential long-term residents, and whether natives have opportunities for meaningful interaction.
The Greek Case and Existing Gaps
Despite its centrality to the refugee crisis, Greece has received little systematic study. In 2015, more than half of all asylum seekers entering Europe passed through Greek islands (UNHCR 2016). Yet the political consequences of this unprecedented inflow remain underexplored. The existing literature leaves two critical gaps.
- Case Gap: Greece, one of the most affected countries, has not been studied in depth, even though it is home to Golden Dawn, the most extreme-right party represented in a European parliament at the time (Heinö 2016).
 - Mechanism Gap: Prior work often conflates exposure and contact. Yet in contexts where refugees are only briefly present—such as the Aegean islands, where most continued to the mainland within 48 hours (Capon 2015)—contact theory cannot apply. This makes it possible to isolate the effect of exposure alone, independent of sustained interaction or direct competition for resources.
 
Contribution
This paper addresses these gaps by studying Greece during the 2015 refugee crisis. It introduces new evidence from a natural experiment in the Aegean, where islands were differentially exposed to refugee arrivals depending on their proximity to Turkey. By using a difference-in-differences strategy (Bertrand, Duflo, and Mullainathan 2004), the study provides causal estimates of refugee exposure on extreme-right support. The contribution is threefold. First, it shows that mere exposure—even without contact or competition—can increase far-right vote share. Second, it clarifies why previous studies have produced contradictory findings: the mechanism of exposure (contact vs. transient presence) matters. Third, it expands the comparative scope by analyzing Greece, a frontline country and unique test case for the “flash potential” of refugee crises to reshape politics (Sniderman, Hagendoorn, and Prior 2004; Hopkins 2010).
3 Theory and Argument
Key Concepts
The dependent variable in this study is electoral support for Golden Dawn (GD), Greece’s extreme-right party that combined ultranationalist rhetoric with anti-immigrant and xenophobic appeals (Heinö 2016). The independent variable of interest is exposure to refugees, defined as the sudden and localized presence of asylum seekers on Greek islands during the 2015 refugee crisis. Importantly, “exposure” here refers to the physical arrival and temporary presence of refugees, rather than long-term settlement or integration.
Theoretical Framework
Two dominant theories in the literature offer competing expectations about how refugee inflows affect political behavior.
- Contact theory argues that intergroup contact under appropriate conditions reduces prejudice (Allport 1979). When natives and refugees interact meaningfully, hostility may decrease.
 - Realistic group conflict theory suggests that competition over scarce resources—jobs, housing, welfare—fuels intergroup hostility and boosts support for exclusionary politics (Campbell 1965).
 
Yet the Greek case challenges the assumptions of both frameworks. Refugees were present on the islands only briefly, often less than 48 hours (Capon 2015). This limited duration precluded sustained interaction (undermining the scope of contact theory) and prevented meaningful resource competition (limiting the applicability of conflict theory).
Argument
This article argues that exposure alone—even without contact or competition—can increase support for the extreme right. The mechanism lies in the symbolic and perceptual shock created by the sudden arrival of large numbers of refugees. Even temporary visibility may be sufficient to trigger perceptions of cultural threat, insecurity, or loss of control. Refugee arrivals can thus activate latent nationalist and xenophobic predispositions, making extremist parties more attractive.
Hypotheses
H1 (Exposure Hypothesis): Municipalities exposed to refugee arrivals will experience an increase in Golden Dawn’s vote share relative to non-exposed municipalities.
H2 (Dose-Response Hypothesis): The effect will be stronger in municipalities with higher per-capita refugee arrivals.
Causal Framework (DAG)
The causal logic can be represented as follows:
- Treatment (IV): Refugee Arrivals (driven by proximity to Turkey)
 - Mediator: Perceptions of Threat (cultural, security, symbolic)
 - Outcome (DV): Golden Dawn Vote Share
 
Show the code
# ------------------------------------------------------------
# Goal: Draw a simple Directed Acyclic Graph (DAG) that represents
# a theory of how exposure to refugees affects Golden Dawn (GD) vote share.
# This DAG uses ggdag + dagitty + ggplot2.
#
# Key idea: We draw *nodes* (circles) for variables and *arrows* for
#           hypothesized causal effects. We'll also show how to make
#           arrows shorter so they don't touch the nodes.
# ------------------------------------------------------------
# --- 0) Playground: tweak these safely -------------------------------
node_size    <- 25    # how big each node (circle) is
node_alpha   <- 0.5   # transparency of node fill (0 = invisible, 1 = solid)
arrowhead_pt <- 10    # arrowhead size in points (this is the *tip* triangle)
trim         <- 0.20  # how much of the arrow *shaft* to chop off from BOTH ends
                     # e.g., 0.20 = remove 20% at the start AND 20% at the end
                     # so the arrow in the middle is 60% of the original length
# --- 1) Load required packages ----------------------------------------------
# ggdag   → helper functions to draw DAGs in ggplot2
# dagitty → defines DAG objects and relationships between variables
# ggplot2 → plotting system used to visualize the DAG
library(ggdag)
library(dagitty)
library(ggplot2)
# --- 2) Define the DAG structure --------------------------------------------
# dagify() creates a DAG object where you specify arrows (causal paths).
# The notation A ~ B means "A is caused by B".
# Here:
#   GDvote   ← threat     (Golden Dawn vote share caused by perceived threat)
#   threat   ← exposure   (Threat is caused by exposure to refugees)
#   exposure ← distance   (Exposure caused by proximity to Turkey)
# Labels: provide readable names for the DAG nodes (with line breaks).
# Coords: x and y positions of each node (for nice plotting layout).
theory_dag <- dagify(
  GDvote ~ threat,             # DV (Golden Dawn vote share) depends on threat
  threat ~ exposure,           # Threat depends on exposure
  exposure ~ distance,         # Exposure depends on distance
  labels = c(
    GDvote   = "Golden Dawn\nVote Share (DV)",
    threat   = "Perceived Threat\n(Mediator)",
    exposure = "Refugee Arrivals",
    distance = "Proximity to Turkey"
  ),
  coords = list(
    x = c(distance = 0, exposure = 1, threat = 2, GDvote = 3),
    y = c(distance = 1, exposure = 1, threat = 1, GDvote = 1)
  )
)
# --- 3) Convert DAG into a data frame for ggplot ----------------------------
# tidy_dagitty() extracts coordinates, labels, and edges for plotting.
# as.data.frame() makes it into a tibble/data frame usable by ggplot.
dag_df <- as.data.frame(tidy_dagitty(theory_dag, layout = "auto"))
# Extract plotting limits (so the nodes don’t touch the plot borders).
min_x <- min(dag_df$x); max_x <- max(dag_df$x)
min_y <- min(dag_df$y); max_y <- max(dag_df$y)
error <- 0.2  # small padding around the plot
# --- 3.5) Make arrow *segments* shorter and floating between nodes ----------
# Important: In ggplot, the "arrow()" only controls the arrowhead (the little
# triangle tip). It does NOT control the length of the shaft. Shaft length is
# determined by the start (x, y) and end (xend, yend) coordinates.
#
# To make arrow shafts shorter AND not start at the node centers, we:
# 1) Find all edge rows (they have xend/yend values).
# 2) Move the start point forward by 'trim' proportion along the edge.
# 3) Move the end point backward by 'trim' proportion along the edge.
# The result: the arrow appears between nodes without touching them.
edges_df <- subset(dag_df, !is.na(xend) & !is.na(yend))
# Safety checks: ensure trim is between 0 and 0.49
# (If trim >= 0.5 the start could pass the end!)
trim <- max(min(trim, 0.49), 0.00)
# Compute trimmed start (xstart_trim, ystart_trim) and end (xend_trim, yend_trim)
edges_df$xstart_trim <- edges_df$x   + trim * (edges_df$xend - edges_df$x)
edges_df$ystart_trim <- edges_df$y   + trim * (edges_df$yend - edges_df$y)
edges_df$xend_trim   <- edges_df$xend - trim * (edges_df$xend - edges_df$x)
edges_df$yend_trim   <- edges_df$yend - trim * (edges_df$yend - edges_df$y)
# --- 4) Plot the DAG --------------------------------------------------------
# Each node is a big semi-transparent point with a label on top.
# Edges (arrows) are drawn with geom_dag_edges_arc using the *trimmed* endpoints.
ggplot(data = dag_df) +
  # Draw big transparent points for each node
  geom_dag_point(aes(x = x, y = y), size = node_size, alpha = node_alpha) +
  
  # Add text labels (only once per node → !duplicated(label))
  geom_label(
    data = subset(dag_df, !duplicated(label)),
    aes(x = x, y = y, label = label),
    fill = alpha("white", 0.8)
  ) +
  
  # Draw arrows between nodes (edges of the DAG)
  # NOTE: arrow(length = ...) controls ONLY the arrowhead size (the tip).
  # The shaft length is controlled by the start/end coordinates we just trimmed.
  geom_dag_edges_arc(
    data = edges_df,
    aes(x = xstart_trim, y = ystart_trim, xend = xend_trim, yend = yend_trim),
    curvature = 0.0,
    arrow = grid::arrow(length = grid::unit(arrowhead_pt, "pt"), type = "closed")
  ) +
  
  # Set x- and y-limits to keep everything inside the plot
  coord_sf(
    xlim = c(min_x - error, max_x + error),
    ylim = c(min_y - error, max_y + error)
  ) +
  
  # Minimal background → removes axes, gridlines, etc.
  theme_void()Scope Conditions
The argument applies under specific conditions:
- Temporal scope: Situations of sudden refugee crises, rather than gradual, long-term migration.
 - Spatial scope: Frontline states where arrivals are geographically concentrated and highly visible (e.g., Greek islands, Italian coastal towns).
 - Political scope: Contexts with viable extreme-right parties able to mobilize around issues of national identity and security.
 
4 Methods
Research Design
This study uses a quantitative research design to estimate the causal effect of refugee arrivals on electoral support for Golden Dawn (GD). The research design exploits the fact that geography determined which Greek islands were exposed to sudden inflows of refugees in 2015. Islands closer to the Turkish coast received large numbers of arrivals, while other islands with similar social and economic characteristics received few or none. This quasi-random exposure provides a natural experiment that allows for causal inference.
The central strategy is a difference-in-differences (DID) design, complemented by an event-study analysis. DID compares changes in GD vote share between exposed and non-exposed municipalities before and after the refugee crisis. The event study further assesses the validity of the parallel trends assumption and the dynamics of the effect over time.
The unit of analysis is the municipality (N = 95). The temporal domain covers four national parliamentary elections between 2012 and 2015, which are coded as four consecutive periods:
- May 2012 = 2012
 - June 2012 = 2013
 - January 2015 = 2015
 - September 2015 = 2016
 
Data
Two primary sources are used:
Electoral outcomes: Official election returns from the Greek Ministry of Interior and Public Administration, at a municipality level, for all four elections.
Refugee arrivals: Monthly data from the UNHCR, aggregated to the municipal level and expressed as arrivals per capita. I truncate the measure at five arrivals per resident (the 99.7th percentile of the non-zero distribution) to prevent extreme values from skewing results.
The dependent variable is the vote share of Golden Dawn at the local level. The independent variable is exposure to refugee arrivals, coded either as a binary treatment (any arrivals vs. none) or as the per capita number of arrivals.
Difference-in-Differences Model
The baseline DID model is specified as:
\[ GD_{st} = \gamma_{s} + \lambda_{t} + \delta_{DID} T_{st} + u_{st} \]
where:
- \(GD_{st}\) is the vote share for Golden Dawn in municipality \(s\) at election \(t\)
 - \(\gamma_{s}\) is a unit fixed effect controlling for time-invariant characteristics of each municipality
 - \(\lambda_{t}\) is an election fixed effect controlling for shocks common to all units in a given election
 - \(T_{st}\) is the treatment indicator (binary or continuous exposure to refugee arrivals)
 - \(u_{st}\) is an idiosyncratic error term
 
The coefficient of interest, \(\delta_{DID}\), captures the causal effect of refugee arrivals on GD vote share.
Event-Study Model
To probe the validity of the DID design and assess dynamic effects, I estimate an event-study specification:
\[ GD_{st} = \gamma_{s} + \lambda_{t} + \sum_{k \neq -1} \delta_{k} \, D_{t=k} \times Treatment_{s} + \alpha_{s} t + u_{st} \]
Here, \(D_{t=k}\) are dummies for election periods, with \(k=-1\) (January 2015) as the baseline. The interaction terms capture how treated and untreated municipalities evolve across time. Unit-specific linear trends (\(\alpha_s t\)) allow each municipality to follow its own trajectory, thereby accounting for unobserved local factors. Event-study plots (see Figure 3) show that treated and untreated areas followed parallel trends before September 2015 and diverged only after the refugee crisis, supporting a causal interpretation.
5 Findings
Exposure and Electoral Support for Golden Dawn
The analysis shows that exposure to refugees significantly increased electoral support for Golden Dawn (GD). Figure 3 plots average GD vote shares in treated and untreated municipalities across the four elections (2012–2016). Both groups followed parallel trends before the refugee crisis, diverging only in September 2015. In that election, municipalities exposed to refugee arrivals saw GD support increase by nearly two percentage points, while support in non-exposed areas remained flat. This provides strong evidence that the observed effect is driven by refugee inflows rather than by pre-existing differences in political trajectories.
Show the code
# ------------------------------------------------------------
# Goal: Prepare panel features (trends, lags, diffs, dummies),
# build grouped summaries with CIs, and plot GD vote % over time
# by treatment status.
# ------------------------------------------------------------
# --- 0) Make sure 'year' is numeric -----------------------------------------
# Many imported datasets store years as character or labelled types.
# Converting to integer avoids problems when sorting, lagging, etc.
datamuni$year <- as.integer(datamuni$year)
# --- 1) Municipality-specific linear trends (Model 2, etc.) -----------------
# Some specifications include a separate linear trend for each municipality.
# Your model expects columns named s2trend, s3trend, ..., s94trend explicitly.
# Each trend is "year" for that municipality, and 0 for all others.
# NOTE:
# - This starts at m = 2 to match your model’s expected naming.
# - Requires 'datamuni$muni' to be an integer ID in [1..94].
for (m in 2:94) {
  trend_col <- paste0("s", m, "trend")
  datamuni[[trend_col]] <- ifelse(datamuni$muni == m, datamuni$year, 0L)
}
# --- 2) Lags (only for GD in percentage terms) -------------------------------
# We create:
#   - gdper     : Golden Dawn vote share in percent
#   - gdlagper  : 1-period lag of GD vote share (in percent)
#   - gdlag3per : 3-period lag (in percent)
# Ordering by muni, then year, ensures lags are computed in time order.
# IMPORTANT:
# - Lags produce NA in the first available years per municipality.
# - 'gdvote' is assumed to be a proportion (0–1). Multiply by 100 for %.
library(dplyr)
datamuni <- datamuni %>%
  arrange(muni, year) %>%
  group_by(muni) %>%
  mutate(
    gdper     = gdvote * 100,                    # main DV in percent
    gdlagper  = dplyr::lag(gdvote,  1) * 100,    # 1-year lag, in percent
    gdlag3per = dplyr::lag(gdvote,  3) * 100     # 3-year lag, in percent
  ) %>%
  ungroup()
# --- 3) Differences for Reduced Form / Placebo figures -----------------------
# First difference of GD % (current minus 1-year lag), used in fig4a.
# A "placebo" style difference using 1-year-lag minus 3-year-lag, used in fig4b.
# NOTE: These will be NA where required lags are missing.
datamuni <- datamuni %>%
  mutate(
    gdperdif  = gdper    - gdlagper,   # ΔGD% (t − t−1)
    gdperdif3 = gdlagper - gdlag3per   # ΔGD% (t−1 − t−3)
  )
# --- 4) Year dummies + interactions (Table S5 / Fig. S4) ---------------------
# Create indicator (dummy) variables for specific years and interact with logdist.
# Here, y1..y4 correspond to years 2012, 2013, 2015, 2016 respectively.
# The interactions (y*dist) let the effect of distance vary by those years.
# Assumes 'logdist' already exists (e.g., log(distance to Turkey)).
datamuni <- datamuni %>%
  mutate(
    y1 = as.integer(year == 2012),
    y2 = as.integer(year == 2013),
    y3 = as.integer(year == 2015),
    y4 = as.integer(year == 2016),
    y1dist = y1 * logdist,
    y2dist = y2 * logdist,
    y3dist = y3 * logdist,
    y4dist = y4 * logdist
  )
# --- 5) Other percent levels for appendix summaries --------------------------
# Convert other commonly used variables to percent levels (no lags/diffs here).
# 'ndvote' and 'turnout' are assumed to be proportions (0–1).
datamuni <- datamuni %>%
  mutate(
    ndper   = ndvote  * 100,   # New Democracy in %
    turnper = turnout * 100    # Turnout in %
  )
# --- 6) Grouped means, SDs, and 95% CIs by treatment & year ------------------
# We compute:
#   - mean GD% (gdper_mean)
#   - SD of GD% (gdper_sd)
#   - sample size (count_n)
#   - standard error (se = sd/sqrt(n))
#   - 95% CI ≈ mean ± 1.96 * se (normal approx)
# NOTE:
# - If your data contain NA values, consider mean(..., na.rm = TRUE) / sd(..., na.rm = TRUE).
# - 'treatment_group' is assumed to be coded (e.g., 0 = Control, 1 = Treated).
df_muni_avg2 <- datamuni %>%
  mutate(gdper = gdvote * 100) %>%
  group_by(treatment_group, year) %>%
  summarize(
    gdper_mean = mean(gdper),   # add na.rm = TRUE if needed
    gdper_sd   = sd(gdper),     # add na.rm = TRUE if needed
    count_n    = n(),
    .groups    = "drop"
  ) %>%
  mutate(
    se_gdper = gdper_sd / sqrt(count_n),
    gdcilo   = gdper_mean - 1.96 * se_gdper,   # lower 95% bound
    gdcihi   = gdper_mean + 1.96 * se_gdper    # upper 95% bound
  )
# --- 7) Plot: GD% over time by treatment group with 95% CIs ------------------
# - Lines/points show average GD% over time within each group.
# - Error bars visualize 95% confidence intervals.
# - Colors: blue = Control, red = Treated (per your manual scale).
# - Y-axis fixed to [2, 11] with 1-pt ticks for readability (match your figure).
fig2a <- ggplot(
  df_muni_avg2,
  aes(x = year, y = gdper_mean,
      color = as.factor(treatment_group),
      group = treatment_group)
) +
  geom_line(size = 1, position = position_dodge(width = 0.2)) +
  geom_point(size = 2.5, position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = gdcilo, ymax = gdcihi),
                width = 0.2,
                position = position_dodge(width = 0.2)) +
  scale_color_manual(
    values = c("blue", "red"),
    labels = c("Control", "Treated"),
    name   = "Group"
  ) +
  theme_bw() +
  scale_y_continuous(
    limits = c(2, 11),
    breaks = seq(2, 11, by = 1)
  ) +
  theme(
    legend.position      = c(0.00, 1.00),   # place legend at top-left inside plot
    legend.justification = c("left", "top"),
    legend.background    = element_rect(fill = "white", color = "black")
  )
# --- (Optional) Save the figure to disk --------------------------------------
# ggsave(here::here("figures", "gd_trends_by_treatment.png"),
#        fig2a, width = 7, height = 4.5, dpi = 300)
# Display the figure
print(fig2a)Main Difference-in-Differences Estimates
Table 1 presents the regression results. Model 1 shows that exposed municipalities experienced a 2 percentage point increase in GD vote share compared to controls. This effect is substantial, representing roughly a 44 percent increase relative to GD’s baseline support (4.5 percent in January 2015). Model 2 confirms that the effect is robust to the inclusion of municipality-specific linear trends, ruling out the possibility that smooth local dynamics explain the findings.
Placebo Tests
To assess whether the effect might be spurious, placebo models use lagged GD vote share as the dependent variable. Models 3 and 4 in Table 1 show no significant association between refugee arrivals and GD support in previous elections, providing reassurance that the September 2015 result is not driven by pre-treatment dynamics.
Intensity of Treatment
Replacing the binary treatment indicator with the per capita number of refugee arrivals yields consistent results. Models 9 and 10 show that a higher ratio of refugees to residents is associated with larger increases in GD support. By contrast, placebo tests using lagged outcomes (Models 11 and 12) show no effect before 2015. This dose-response pattern further supports the causal interpretation.
Show the code
# ------------------------------------------------------------
# Goal: Estimate fixed effects OLS models of Golden Dawn vote share
# using different treatment variables (binary vs. arrivals proportion),
# with/without municipality-specific trends, and report results in a table.
# ------------------------------------------------------------
# --- 1) Load packages --------------------------------------------------------
# fixest → fast and flexible fixed effects regressions
# broom  → convert model output to a tidy data frame
# dplyr  → for filtering results
# modelsummary → for producing nice regression tables
library(fixest)
library(broom)
library(dplyr)
library(modelsummary)
# --- 2) Municipality-level models -------------------------------------------
# Dataset: datamuni (panel of municipalities × years)
# Outcomes: 
#   - gdper     : GD vote share in % (DV)
#   - gdlagper  : lagged GD vote share in %
# Main explanatory variables:
#   - treatment  : binary treatment indicator
#   - trarrprop  : refugee arrivals per capita (continuous treatment)
# FE structure: 
#   - | muni  → municipality fixed effects
# Time dummies: 
#   - i(year) or factor(year)
# SEs: 
#   - cluster = ~muni (robust at municipality level)
###########################
# Models 1, 2, 3, 4, 9, 10, 11, 12
###########################
# Model 1: GD% on binary treatment + year FE + muni FE; cluster by muni
model1 <- feols(gdper ~ i(year) + treatment | muni, 
                data = datamuni, cluster = ~muni)
coef1 <- tidy(model1) %>% filter(term == "treatment")
# Model 2: Add municipality-specific linear trends (s2trend ... s94trend)
trend_vars <- paste0("s", 2:94, "trend")
fmla <- as.formula(
  paste("gdper ~ i(year) + treatment +", 
        paste(trend_vars, collapse = " + "), 
        "| muni")
)
model2 <- feols(fmla, data = datamuni, cluster = ~muni)
coef2 <- tidy(model2) %>% filter(term == "treatment")
# Model 3: Lagged DV (gdlagper) on treatment + FE
model3 <- feols(gdlagper ~ i(year) + treatment | muni,  
                data = datamuni, cluster = ~muni)
coef3 <- tidy(model3) %>% filter(term == "treatment")
# Model 4: Lagged DV + treatment + muni trends
fmla <- as.formula(
  paste("gdlagper ~ i(year) + treatment +", 
        paste(trend_vars, collapse = " + "), 
        "| muni")
)
model4 <- feols(fmla, data = datamuni, cluster = ~muni)
coef4 <- tidy(model4) %>% filter(term == "treatment")
# Model 9: Replace treatment with arrivals per capita (trarrprop)
model9 <- feols(gdper ~ i(year) + trarrprop | muni, 
                data = datamuni, cluster = ~muni)
coef9 <- tidy(model9) %>% filter(term == "trarrprop")
# Model 10: Arrivals per capita + muni trends
fmla2 <- as.formula(
  paste("gdper ~ trarrprop + i(year) +", 
        paste(trend_vars, collapse = " + "), 
        "| muni")
)
model10 <- feols(fmla2, data = datamuni, cluster = ~muni)
coef10 <- tidy(model10) %>% filter(term == "trarrprop")
# Model 11: Lagged DV (gdlagper) + arrivals per capita
model11 <- feols(gdlagper ~ i(year) + trarrprop | muni, 
                 data = datamuni, cluster = ~muni)
coef11 <- tidy(model11) %>% filter(term == "trarrprop")
# Model 12: Lagged DV + arrivals per capita + muni trends
fmla <- as.formula(
  paste("gdlagper ~ factor(year) + trarrprop +", 
        paste(trend_vars, collapse = " + "), 
        "| muni")
)
model12 <- feols(fmla, data = datamuni, cluster = ~muni)
coef12 <- tidy(model12) %>% filter(term == "trarrprop")
# --- 3) Notes ---------------------------------------------------------------
# - i(year) and factor(year) both add year dummies (fixed effects for time).
# - tidy() + filter(term == "...") pulls out only the coefficient of interest.
# - Clustered SEs should generally match the FE unit (here: municipality).
# - Trends s2trend...s94trend were created earlier.
# --- 4) Produce regression table --------------------------------------------
# Create a named list of models to display side by side
models <- list(
  "(1)" = model1, "(2)" = model2, "(3)" = model3, "(4)" = model4,
  "(9)" = model9, "(10)" = model10, "(11)" = model11, "(12)" = model12
)
# Use modelsummary to display results
modelsummary(
  models,
  coef_map = c(
    treatment = "Binary treatment",
    trarrprop = "Arrivals per capita"
  ),
  estimate   = "{estimate}{stars}",   # show estimate with significance stars
  statistic  = "({std.error})",       # show robust SEs in parentheses
  stars      = TRUE,                  # add stars for p < 0.1, 0.05, 0.01
  gof_omit   = "AIC|BIC|Log|R2",      # omit some goodness-of-fit stats
  output     = "html"                 # can also use "latex", "markdown", etc.
)| (1) | (2) | (3) | (4) | (9) | (10) | (11) | (12) | |
|---|---|---|---|---|---|---|---|---|
| Binary treatment | 2.079*** | 2.112** | −0.040 | −0.055 | ||||
| (0.351) | (0.674) | (0.392) | (0.713) | |||||
| Arrivals per capita | 0.006 | 0.002 | 0.001 | −0.009*** | ||||
| (0.004) | (0.005) | (0.001) | (0.002) | |||||
| Num.Obs. | 380 | 380 | 285 | 285 | 379 | 379 | 284 | 284 | 
| RMSE | 0.95 | 0.61 | 0.95 | 0.44 | 1.00 | 0.64 | 0.95 | 0.43 | 
| Std.Errors | by: muni | by: muni | by: muni | by: muni | by: muni | by: muni | by: muni | by: muni | 
| FE: muni | X | X | X | X | X | X | X | X | 
Robustness: Distance to Turkey
The reduced-form evidence reinforces these findings. Figure 4 shows that the change in GD vote share between January and September 2015 is strongly correlated with proximity to Turkey. By contrast, the placebo comparison (2012–2015) shows no such relationship. Appendix Table S2 and Figure S1 confirm that distance predicted GD support only in September 2015, consistent with the causal role of refugee arrivals.
Show the code
# ------------------------------------------------------------
# Goal: Plot reduced-form and placebo relationships between
# logged distance to Turkey and changes in Golden Dawn vote share.
# ------------------------------------------------------------
# --- 1) Subset the data to 2016 only ----------------------------------------
# We focus on 2016 (the election of interest).
# 'datamuni' must already contain:
#   - logdist     : logged distance to Turkey
#   - gdperdif    : ΔGD% (current − lagged, Jan–Sep 2015 → main reduced form)
#   - gdperdif3   : placebo difference (lagged − 3-lagged)
red <- subset(datamuni, year == 2016)
# --- 2) Reduced Form Figure -------------------------------------------------
# Scatterplot: x = log distance, y = change in GD vote share
# Add linear fit (stat_smooth with method = "lm")
# Limit axes so the figure matches published version (y: −5..5, x: 0.2..6.3)
fig4a <- ggplot(red, aes(x = logdist, y = gdperdif)) + 
  geom_point(shape = 1) +                           # open circles
  ylim(-5, 5) +                                     # y-axis limits
  xlim(0.2, 6.3) +                                  # x-axis limits
  stat_smooth(method = "lm") +                      # add OLS fit line
  xlab("Logged Distance in Km") +                   # x-axis label
  ylab("Change in GD vote: Jan–Sep 2015") +         # y-axis label
  theme(legend.position = "none") +                 # no legend needed
  ggtitle("Reduced Form") +                         # title
  theme_bw()                                        # clean theme
# --- 3) Placebo Figure ------------------------------------------------------
# Same plot setup, but outcome is gdperdif3 (placebo difference).
# The placebo test checks whether “fake” changes line up with distance.
fig4b <- ggplot(red, aes(x = logdist, y = gdperdif3)) + 
  geom_point(shape = 1) +
  ylim(-5, 5) +
  xlim(0.2, 6.3) +
  stat_smooth(method = "lm") +
  xlab("Logged Distance in Km") +
  ylab("Change in GD vote: Jan–Sep 2015") +
  theme(legend.position = "none") +
  ggtitle("Placebo") +
  theme_bw()
# --- 4) Arrange the two figures side by side --------------------------------
# ggarrange() puts both figures into one row with two columns.
ggarrange(fig4a, fig4b, ncol = 2)Additional Robustness
The appendix further strengthens the results. Table S1 reports descriptive statistics, showing comparable pre-crisis distributions between treated and control municipalities. Figure S2 demonstrates that the estimated treatment effect remains stable across a wide range of distance cutoffs, alleviating concerns that islands close to Turkey differ systematically from distant ones.
Summary of Findings
Taken together, these results provide clear support for the core hypothesis: brief exposure to refugees, even without long-term settlement or competition, can increase support for extreme-right parties. The DID estimates establish a sharp causal effect, placebo tests confirm validity, and dose-response as well as geographic proximity analyses reinforce the robustness of the results.
In the broader literature, these findings align with Steinmayr’s (Steinmayr 2016) evidence from Austria that exposure can shape voting, but they diverge from Dustmann et al. (Dustmann, Vasiljeva, and Damm 2016), who found that contact reduced hostility in Denmark. Greece thus represents a crucial case: where exposure is sudden, highly visible, and short-lived, it increases rather than decreases support for the extreme right.
6 Discussion
Contribution
The findings show that brief and localized exposure to refugees increased electoral support for Golden Dawn, Greece’s extreme-right party. This challenges the dominant focus in the literature on long-term labor migration and resource competition (Barone et al. 2016; Halla, Wagner, and Zweimüller 2015; Becker and Fetzer 2016). Instead, the results highlight how even short-lived encounters—without sustained contact or economic rivalry—can shift political outcomes. In doing so, the study contributes to debates about the “flash potential” of refugee crises (Sniderman, Hagendoorn, and Prior 2004; Hopkins 2010) and demonstrates that refugee inflows must be understood as distinct political shocks rather than extensions of general immigration dynamics.
Broader Implications
These results underline the symbolic and affective dimensions of immigration politics. Refugees’ brief passage through Aegean islands was sufficient to trigger electoral backlash, suggesting that perceptions of disorder, security, and national identity may matter as much as economic or demographic pressures. This complicates contact and conflict theories (Allport 1979; Campbell 1965), which assume more sustained interaction. It also points to the need for policymakers and scholars to treat refugee arrivals as unique events with the potential to rapidly reshape the political landscape.
Limitations and Future Research
The analysis has three main limitations. First, the study focuses on Greek islands in 2015, so generalizability to mainland areas or other contexts may be limited. Second, while the quasi-experimental design strengthens causal claims, it does not fully identify the mechanisms at work—whether media coverage, perceived loss of control, or nationalist mobilization. Third, the analysis centers on electoral outcomes and cannot capture broader civic or attitudinal responses. Future research should therefore combine cross-national comparisons, survey evidence, and media analysis to explore mechanisms and scope conditions.
7 Conclusion
Restating the Argument
This article has demonstrated that localized exposure to the 2015 refugee crisis causally increased electoral support for Golden Dawn, Greece’s extreme-right party. Even when contact with refugees was brief and resource competition minimal, municipalities receiving arrivals recorded a significant rise in Golden Dawn’s vote share.
Contribution to the Literature
These findings contribute to debates on immigration and political behavior by showing that refugees function as a distinct type of political shock. Unlike labor migrants, whose long-term presence tends to shape party competition through economic channels, refugees can trigger immediate political reactions grounded in perceptions of threat and disorder. In doing so, the analysis pushes beyond existing studies of Austria and Denmark to highlight the specific dynamics of frontline states.
Broader Implications
For both scholars and a general audience, the findings underscore the fragility of political systems to sudden demographic disruptions. The refugee crisis was not only a humanitarian and economic challenge but also a catalyst for populist and extreme-right mobilization. This perspective reframes how we think about the link between migration crises and democratic resilience.
Unresolved Questions
Important gaps remain. What mechanisms explain why short-lived exposure produced lasting electoral consequences—media framing, nationalist rhetoric, or voter perceptions of state weakness? Under what conditions might such effects fade, rather than entrench? Future research should extend this analysis cross-nationally, linking electoral outcomes with survey and media data to clarify how refugee crises shape political trajectories.
8 References
9 Appendix
9.1 Summary Statistics
Table S1 displays the descriptive statistics of all variables used in the analysis at the municipality level.
Show the code
# ------------------------------------------------------------
# Goal: Produce a summary statistics table for key variables
# across different election years and treatments.
# ------------------------------------------------------------
# --- 1) Load required packages ----------------------------------------------
# dplyr  → data wrangling
# purrr  → functional programming tools (imap_dfr to loop + label)
# tibble → cleaner data frame format
library(dplyr)
library(purrr)
library(tibble)
# --- 2) Create year-specific subsets ----------------------------------------
# Focus on a few key election periods. Keep only rows with valid muni ID.
# Note: year codes follow how the dataset is structured:
#   - 2016 subset = used for Sep 2015 election outcome
#   - 2015 subset = Jan 2015 election
#   - 2013 subset = Jun 2012 election
#   - 2012 subset = Jan 2012 election
sum_df_sep_2015 <- subset(datamuni, year == 2016 & !is.na(muni))
sum_df_jan_2015 <- subset(datamuni, year == 2015 & !is.na(muni))
sum_df_jan_2012 <- subset(datamuni, year == 2012 & !is.na(muni))
sum_df_jun_2012 <- subset(datamuni, year == 2013 & !is.na(muni))
# --- 3) Define variables to summarize ---------------------------------------
# We build a named list: each entry = "Label" -> data vector.
# Labels will appear in the summary table.
spec <- list(
  "Binary treatment"                         = sum_df_sep_2015$treatment,
  "Arrivals per capita"                      = sum_df_sep_2015$trarrprop,
  "Log distance"                             = sum_df_sep_2015$logdist,
  "Turnout (%) (Sep 2015)"                   = sum_df_sep_2015$turnper,
  "Turnout (%) (Jan 2015)"                   = sum_df_jan_2015$turnper,
  "GD vote share (%) (Sep 2015)"             = sum_df_sep_2015$gdper,
  "GD vote share (%) (Jan 2015)"             = sum_df_jan_2015$gdper,
  "GD vote share (%) (Jun 2012)"             = sum_df_jun_2012$gdper,
  "GD vote share (%) (Jan 2012)"             = sum_df_jan_2012$gdper,
  "Nea Dimokratia vote share (%) (Sep 2015)" = sum_df_sep_2015$ndper,
  "Nea Dimokratia vote share (%) (Jan 2015)" = sum_df_jan_2015$ndper
)
# --- 4) Compute summary statistics ------------------------------------------
# imap_dfr loops over (value, name) pairs in the list `spec`.
# For each variable vector (x):
#   - N   = number of non-missing obs
#   - Mean, SD, Min, Max = rounded to 3 decimals
# Results are stacked into a single tibble.
summary_df <- imap_dfr(spec, ~{
  x <- .x  # data vector
  tibble(
    Variable = .y,
    N        = sum(!is.na(x)),
    Mean     = round(mean(x, na.rm = TRUE), 3),
    SD       = round(sd(x,   na.rm = TRUE), 3),
    Min      = round(min(x,  na.rm = TRUE), 3),
    Max      = round(max(x,  na.rm = TRUE), 3)
  )
})
# --- 5) Display as a nice table ---------------------------------------------
# knitr::kable formats the tibble into a readable table (HTML here).
knitr::kable(
  summary_df,
  format = "html",   # use "latex" if compiling to PDF
  caption = "Summary statistics of the dataset"
)| Variable | N | Mean | SD | Min | Max | 
|---|---|---|---|---|---|
| Binary treatment | 95 | 0.126 | 0.334 | 0.000 | 1.000 | 
| Arrivals per capita | 94 | 1.613 | 12.926 | 0.000 | 125.311 | 
| Log distance | 95 | 4.861 | 1.240 | 0.588 | 6.235 | 
| Turnout (%) (Sep 2015) | 95 | 48.405 | 9.757 | 18.035 | 64.035 | 
| Turnout (%) (Jan 2015) | 95 | 51.570 | 13.353 | 17.791 | 74.319 | 
| GD vote share (%) (Sep 2015) | 95 | 6.018 | 2.475 | 0.000 | 17.513 | 
| GD vote share (%) (Jan 2015) | 95 | 4.460 | 2.129 | 0.000 | 11.927 | 
| GD vote share (%) (Jun 2012) | 95 | 5.435 | 2.366 | 0.483 | 12.616 | 
| GD vote share (%) (Jan 2012) | 95 | 5.585 | 3.151 | 0.794 | 18.839 | 
| Nea Dimokratia vote share (%) (Sep 2015) | 95 | 28.513 | 6.620 | 12.318 | 47.500 | 
| Nea Dimokratia vote share (%) (Jan 2015) | 95 | 29.836 | 8.334 | 13.093 | 54.145 | 
9.2 Intention to Treat Test
I examine “Intention to Treat” in Table S2. I specifically interact logged distance with each election-dummy (using the May 2012 election as baseline). I show that distance to the Turkish coast predicts GD vote share only in September 2015, but not in any of the previous elections: May 2012, June 2012 and January 2015 elections.
Show the code
library(fixest)
library(modelsummary)
# ------------------------------------------------------------
# Goal: Estimate a model with interactions between log distance
# and year dummies, then display results in a regression table.
# ------------------------------------------------------------
# --- 1) Run regression ------------------------------------------------------
# feols() syntax: outcome ~ predictors, data = ...
# Outcome:
#   - gdper : Golden Dawn vote share in %
# Predictors:
#   - y2, y3, y4 : year dummies (baseline = 2012 Jan)
#   - logdist    : log distance to Turkey
#   - y2dist, y3dist, y4dist : interactions (log distance × year dummies)
# Note: y2–y4 and interaction columns should already exist in datamuni.
mod <- feols(
  gdper ~ y2 + y3 + y4 + logdist + y2dist + y3dist + y4dist,
  data = datamuni
)
# --- 2) Create a list of models ---------------------------------------------
# Useful if you want to report multiple specifications side by side.
# Here we only have one model, labeled "(1)" in the table.
models <- list("(1)" = mod)
# --- 3) Display results in a formatted table --------------------------------
# modelsummary() nicely formats regression results.
# Key options:
#   - coef_map   : renames predictors to human-readable labels
#   - estimate   : "{estimate}{stars}" prints coef + sig stars
#   - statistic  : shows SEs in parentheses
#   - stars      : significance levels added automatically
#   - gof_omit   : drop some fit stats (AIC, BIC, etc.)
#   - output     : format of table ("html" here; can also do "latex")
modelsummary(
  models,
  coef_map = c(
    "logdist"     = "(Logged) Distance",
    "y2"          = "2012 June",
    "y3"          = "2015 January",
    "y4"          = "2015 September",
    "y2dist"      = "Log Distance × 2012 June",
    "y3dist"      = "Log Distance × 2015 January",
    "y4dist"      = "Log Distance × 2015 September",
    "(Intercept)" = "Constant"
  ),
  estimate  = "{estimate}{stars}",    # coefficients + stars
  statistic = "({std.error})",        # robust SEs in parentheses
  stars     = TRUE,                   # show sig. stars
  gof_omit  = "AIC|BIC|Log|R2|RMSE|Std.Errors", # cleaner output
  output    = "html"                  # output type
)| (1) | |
|---|---|
| (Logged) Distance | −0.204 | 
| (0.210) | |
| 2012 June | −1.074 | 
| (1.492) | |
| 2015 January | −0.893 | 
| (1.492) | |
| 2015 September | 2.769+ | 
| (1.492) | |
| Log Distance × 2012 June | 0.190 | 
| (0.298) | |
| Log Distance × 2015 January | −0.048 | 
| (0.298) | |
| Log Distance × 2015 September | −0.481 | 
| (0.298) | |
| Constant | 6.574*** | 
| (1.055) | |
| Num.Obs. | 380 | 
9.3 Intention to Treat Figure
Figure S1 examines “Intention to Treat” displayed in Table S2. Instead of looking directly at how many refugees arrived, I examine the distance to the Turkish coast — because distance is what determined whether an island was likely to receive refugees in the first place.
To do this, I interact (combine) the distance variable with each election period. Think of it as asking:
- In May 2012, was distance related to Golden Dawn’s vote share?
 - In June 2012, was distance related?
 - In January 2015, was distance related?
 - And finally, in September 2015, was distance related?
 
If distance mattered even before refugees arrived, that would be suspicious: maybe islands closer to Turkey were already more right-wing. But the results show the opposite:
- In 2012 and January 2015, distance didn’t predict Golden Dawn’s vote share at all. Islands near and far from Turkey looked similar in terms of voting.
 - In September 2015, right after the refugee crisis, distance suddenly became predictive: closer islands (that actually received the refugees) gave more votes to Golden Dawn.
 
That’s exactly what you’d expect if the effect comes from the refugee arrivals, not from some pre-existing difference. So, Figure S1 shows that distance to the Turkish coast predicts GD vote share only in September 2015, but not in any of the previous elections: May 2012, June 2012 and January 2015 elections.
Show the code
# ------------------------------------------------------------
# Goal: Estimate a regression with distance × year interactions,
# test linear combinations of coefficients (e.g., marginal effect
# of distance in each year), and plot results with 95% CIs.
# ------------------------------------------------------------
# --- 1) Load required packages ----------------------------------------------
library(lmtest)     # for coeftest() and hypothesis tests
library(sandwich)   # for cluster-robust variance-covariance estimators
library(car)        # for linearHypothesis() (testing linear combos)
library(dplyr)      # for mutate() and bind_rows()
library(broom)      # for tidy() to convert results to data frames
library(ggplot2)    # for plotting
# --- 2) Fit the regression model --------------------------------------------
# Outcome:
#   gdper : Golden Dawn vote share (%)
# Predictors:
#   y2, y3, y4 : year dummies (baseline = 2012 May)
#   logdist    : log distance to Turkey
#   y2dist...y4dist : interaction terms (distance × year dummies)
mod <- lm(gdper ~ y2 + y3 + y4 + logdist + y2dist + y3dist + y4dist,
          data = datamuni)
# --- 3) Compute cluster-robust SEs ------------------------------------------
# Cluster by municipality, so SEs allow for within-muni correlation.
vc <- vcovCL(mod, cluster = ~ muni)
# --- 4) Test linear combinations of coefficients ----------------------------
# Idea: In interacted models, the marginal effect of logdist
# differs by year. For example:
#   Effect in 2015 Sep = logdist + y4dist
# We test whether those effects are zero.
## Test: Effect of logdist in 2015 September
lc1 <- linearHypothesis(mod, "logdist + y4dist = 0", vcov. = vc)
lc1_df <- tidy(lc1) %>%
  slice(1) %>%   # keep first row (F-test collapsed into one row)
  mutate(
    min95 = estimate - 1.96 * std.error,  # lower 95% CI
    max95 = estimate + 1.96 * std.error,  # upper 95% CI
    id1   = 1                             # identifier
  )
## Test: Effect of logdist in 2015 January
lc2 <- linearHypothesis(mod, "logdist + y3dist = 0", vcov. = vc)
lc2_df <- tidy(lc2) %>%
  slice(1) %>%
  mutate(
    min95 = estimate - 1.96 * std.error,
    max95 = estimate + 1.96 * std.error,
    id1   = 2
  )
## Test: Effect of logdist in 2012 June
lc3 <- linearHypothesis(mod, "logdist + y2dist = 0", vcov. = vc)
lc3_df <- tidy(lc3) %>%
  slice(1) %>%
  mutate(
    min95 = estimate - 1.96 * std.error,
    max95 = estimate + 1.96 * std.error,
    id1   = 3
  )
## Test: Effect of logdist in 2012 May (baseline year, no interaction)
lc4 <- linearHypothesis(mod, "logdist = 0", vcov. = vc)
lc4_df <- tidy(lc4) %>%
  slice(1) %>%
  mutate(
    min95 = estimate - 1.96 * std.error,
    max95 = estimate + 1.96 * std.error,
    id1   = 4
  )
# --- 5) Combine results into one dataset ------------------------------------
itt <- bind_rows(lc1_df, lc2_df, lc3_df, lc4_df)
# --- 6) Add extra IDs / group labels for plotting ---------------------------
itt <- itt %>%
  mutate(
    id3   = ifelse(id1 %in% 2:4, 2, id1),   # collapse into 2 groups (red vs blue)
    id2   = 5 - id1,                        # reverse order for plotting (y-axis)
    group = ifelse(id3 == 1, "red", "blue") # assign colors
  )
# Convert to factor for ggplot
itt$group <- as.factor(itt$id3)
# --- 7) Plot results --------------------------------------------------------
# Plot estimates with 95% CI error bars.
# - Horizontal error bars show confidence intervals.
# - Vertical dashed line at 0 shows “no effect”.
# - Y-axis labels = election waves.
figs4 <- ggplot(itt, aes(x = estimate, y = id2)) +
  geom_errorbarh(aes(xmax = max95, xmin = min95), height = .1) +
  geom_point(aes(color = group), size = 2.5) +
  scale_color_manual(values = c("red", "blue"), guide = "none") +
  scale_x_continuous(limits = c(-2, 2)) +
  geom_vline(xintercept = 0, linetype = 2, color = "black") +
  xlab("Treatment Effect") + ylab("") +
  scale_y_continuous(
    breaks = c(1, 2, 3, 4),
    labels = c("2012 May", "2012 June", "2015 January", "2015 September")
  ) +
  theme_bw()
# Display figure
figs49.4 Proximity to Turkey Analysis
It is problematic to assume that islands very close to the Turkish coast are similar to islands far from the Turkish coast. Distance ranges from a few kilometers to 530 kilometers.
Figure S2 shows how the effect of refugee exposure remains the same irrespective of whether one focuses on islands in close proximity to the border versus far. I use a DID as in the main analysis.
Fig S7 Data Prep: Show the code
# ------------------------------------------------------------
# Goal: Check robustness to different maximum distances from
# the Turkish coast. For each cutoff (50, 100, ..., 500 km):
#   1) Filter municipalities within the cutoff
#   2) Run DID with muni FE (with and without s-trends)
#   3) Extract the treatment coefficient + 95% CI
#   4) Combine results and plot estimate vs. cutoff
# ------------------------------------------------------------
# --- 1) Load packages --------------------------------------------------------
library(fixest)   # feols(): fixed-effects regressions
library(dplyr)    # filter(), mutate(), bind_rows()
library(broom)    # tidy(): turn model output into data frames
library(ggplot2)  # plotting
# --- 2) Prep: district alias (kept as in your code) --------------------------
# Not used below but retained for consistency with upstream code.
datamuni$district <- datamuni$perfecture
# --- 3) Pre-built municipality trend variable names --------------------------
# Your models include explicit linear trends s2trend ... s94trend.
# We collect those column names for use in formulas below.
trend_vars <- c(
  "s2trend","s3trend","s4trend","s5trend","s6trend","s7trend","s8trend","s9trend",
  "s10trend","s11trend","s12trend","s13trend","s14trend","s15trend","s16trend",
  "s17trend","s18trend","s19trend","s20trend","s21trend","s22trend","s23trend",
  "s24trend","s25trend","s26trend","s27trend","s28trend","s29trend","s30trend",
  "s31trend","s32trend","s33trend","s34trend","s35trend","s36trend","s37trend",
  "s38trend","s39trend","s40trend","s41trend","s42trend","s43trend","s44trend",
  "s45trend","s46trend","s47trend","s48trend","s49trend","s50trend","s51trend",
  "s52trend","s53trend","s54trend","s55trend","s56trend","s57trend","s58trend",
  "s59trend","s60trend","s61trend","s62trend","s63trend","s64trend","s65trend",
  "s66trend","s67trend","s68trend","s69trend","s70trend","s71trend","s72trend",
  "s73trend","s74trend","s75trend","s76trend","s77trend","s78trend","s79trend",
  "s80trend","s81trend","s82trend","s83trend","s84trend","s85trend","s86trend",
  "s87trend","s88trend","s89trend","s90trend","s91trend","s92trend","s93trend","s94trend"
)
# --- 4) Create an empty results data frame -----------------------------------
# We'll append rows (two per cutoff: DD and DD + s-trends).
dist_results <- data.frame()
# ---------------------------- CUT = 50 km ------------------------------------
# Filter to municipalities within 50 km
df50 <- dplyr::filter(datamuni, distance < 50)
# (a) DID with muni FE (no s-trends)
mod_dd_50 <- feols(gdper ~ i(year) + treatment | muni, data = df50, cluster = ~muni)
res_dd_50 <- tidy(mod_dd_50, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 50, group = "DD Estimates")
# (b) DID with muni FE + explicit muni trends
mod_ddstr_50 <- feols(
  as.formula(paste("gdper ~ i(year) + treatment +", paste(trend_vars, collapse = " + "), "| muni")),
  data = df50, cluster = ~muni
)
res_ddstr_50 <- tidy(mod_ddstr_50, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 50, group = "DD with s-trends\nEstimates")
# Add to master results
dist_results <- bind_rows(dist_results, res_dd_50, res_ddstr_50)
# ---------------------------- CUT = 100 km -----------------------------------
df100 <- dplyr::filter(datamuni, distance < 100)
mod_dd_100 <- feols(gdper ~ i(year) + treatment | muni, data = df100, cluster = ~muni)
res_dd_100 <- tidy(mod_dd_100, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 100, group = "DD Estimates")
mod_ddstr_100 <- feols(
  as.formula(paste("gdper ~ i(year) + treatment +", paste(trend_vars, collapse = " + "), "| muni")),
  data = df100, cluster = ~muni
)
res_ddstr_100 <- tidy(mod_ddstr_100, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 100, group = "DD with s-trends\nEstimates")
dist_results <- bind_rows(dist_results, res_dd_100, res_ddstr_100)
# ---------------------------- CUT = 150 km -----------------------------------
df150 <- dplyr::filter(datamuni, distance < 150)
mod_dd_150 <- feols(gdper ~ i(year) + treatment | muni, data = df150, cluster = ~muni)
res_dd_150 <- tidy(mod_dd_150, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 150, group = "DD Estimates")
mod_ddstr_150 <- feols(
  as.formula(paste("gdper ~ i(year) + treatment +", paste(trend_vars, collapse = " + "), "| muni")),
  data = df150, cluster = ~muni
)
res_ddstr_150 <- tidy(mod_ddstr_150, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 150, group = "DD with s-trends\nEstimates")
dist_results <- bind_rows(dist_results, res_dd_150, res_ddstr_150)
# ---------------------------- CUT = 200 km -----------------------------------
df200 <- dplyr::filter(datamuni, distance < 200)
mod_dd_200 <- feols(gdper ~ i(year) + treatment | muni, data = df200, cluster = ~muni)
res_dd_200 <- tidy(mod_dd_200, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 200, group = "DD Estimates")
mod_ddstr_200 <- feols(
  as.formula(paste("gdper ~ i(year) + treatment +", paste(trend_vars, collapse = " + "), "| muni")),
  data = df200, cluster = ~muni
)
res_ddstr_200 <- tidy(mod_ddstr_200, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 200, group = "DD with s-trends\nEstimates")
dist_results <- bind_rows(dist_results, res_dd_200, res_ddstr_200)
# ---------------------------- CUT = 250 km -----------------------------------
df250 <- dplyr::filter(datamuni, distance < 250)
mod_dd_250 <- feols(gdper ~ i(year) + treatment | muni, data = df250, cluster = ~muni)
res_dd_250 <- tidy(mod_dd_250, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 250, group = "DD Estimates")
mod_ddstr_250 <- feols(
  as.formula(paste("gdper ~ i(year) + treatment +", paste(trend_vars, collapse = " + "), "| muni")),
  data = df250, cluster = ~muni
)
res_ddstr_250 <- tidy(mod_ddstr_250, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 250, group = "DD with s-trends\nEstimates")
dist_results <- bind_rows(dist_results, res_dd_250, res_ddstr_250)
# ---------------------------- CUT = 300 km -----------------------------------
df300 <- dplyr::filter(datamuni, distance < 300)
mod_dd_300 <- feols(gdper ~ i(year) + treatment | muni, data = df300, cluster = ~muni)
res_dd_300 <- tidy(mod_dd_300, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 300, group = "DD Estimates")
mod_ddstr_300 <- feols(
  as.formula(paste("gdper ~ i(year) + treatment +", paste(trend_vars, collapse = " + "), "| muni")),
  data = df300, cluster = ~muni
)
res_ddstr_300 <- tidy(mod_ddstr_300, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 300, group = "DD with s-trends\nEstimates")
dist_results <- bind_rows(dist_results, res_dd_300, res_ddstr_300)
# ---------------------------- CUT = 350 km -----------------------------------
df350 <- dplyr::filter(datamuni, distance < 350)
mod_dd_350 <- feols(gdper ~ i(year) + treatment | muni, data = df350, cluster = ~muni)
res_dd_350 <- tidy(mod_dd_350, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 350, group = "DD Estimates")
mod_ddstr_350 <- feols(
  as.formula(paste("gdper ~ i(year) + treatment +", paste(trend_vars, collapse = " + "), "| muni")),
  data = df350, cluster = ~muni
)
res_ddstr_350 <- tidy(mod_ddstr_350, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 350, group = "DD with s-trends\nEstimates")
dist_results <- bind_rows(dist_results, res_dd_350, res_ddstr_350)
# ---------------------------- CUT = 400 km -----------------------------------
df400 <- dplyr::filter(datamuni, distance < 400)
mod_dd_400 <- feols(gdper ~ i(year) + treatment | muni, data = df400, cluster = ~muni)
res_dd_400 <- tidy(mod_dd_400, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 400, group = "DD Estimates")
mod_ddstr_400 <- feols(
  as.formula(paste("gdper ~ i(year) + treatment +", paste(trend_vars, collapse = " + "), "| muni")),
  data = df400, cluster = ~muni
)
res_ddstr_400 <- tidy(mod_ddstr_400, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 400, group = "DD with s-trends\nEstimates")
dist_results <- bind_rows(dist_results, res_dd_400, res_ddstr_400)
# ---------------------------- CUT = 450 km -----------------------------------
df450 <- dplyr::filter(datamuni, distance < 450)
mod_dd_450 <- feols(gdper ~ i(year) + treatment | muni, data = df450, cluster = ~muni)
res_dd_450 <- tidy(mod_dd_450, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 450, group = "DD Estimates")
mod_ddstr_450 <- feols(
  as.formula(paste("gdper ~ i(year) + treatment +", paste(trend_vars, collapse = " + "), "| muni")),
  data = df450, cluster = ~muni
)
res_ddstr_450 <- tidy(mod_ddstr_450, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 450, group = "DD with s-trends\nEstimates")
dist_results <- bind_rows(dist_results, res_dd_450, res_ddstr_450)
# ---------------------------- CUT = 500 km -----------------------------------
df500 <- dplyr::filter(datamuni, distance < 500)
mod_dd_500 <- feols(gdper ~ i(year) + treatment | muni, data = df500, cluster = ~muni)
res_dd_500 <- tidy(mod_dd_500, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 500, group = "DD Estimates")
mod_ddstr_500 <- feols(
  as.formula(paste("gdper ~ i(year) + treatment +", paste(trend_vars, collapse = " + "), "| muni")),
  data = df500, cluster = ~muni
)
res_ddstr_500 <- tidy(mod_ddstr_500, conf.int = TRUE) %>%
  dplyr::filter(term == "treatment") %>%
  dplyr::mutate(distance = 500, group = "DD with s-trends\nEstimates")
dist_results <- bind_rows(dist_results, res_dd_500, res_ddstr_500)
# --- 5) Plot estimates vs. cutoff -------------------------------------------
# Points = treatment effect estimates. Error bars = 95% CI.
# Color distinguishes the two specifications (with vs. without s-trends).
figs7 <- ggplot(dist_results, aes(x = distance, y = estimate, color = group)) +
  geom_point(size = 2, position = position_dodge(24)) +
  scale_color_manual(
    name   = "Estimation\nMethod",
    values = c("red", "blue")
  ) +
  geom_errorbar(
    aes(ymax = conf.high, ymin = conf.low),
    width = 0.3,
    position = position_dodge(width = 24)
  ) +
  scale_x_continuous(
    name   = "Maximum Distance from Turkish Coast in Km.",
    limits = c(40, 520),
    breaks = seq(50, 600, 50)
  ) +
  scale_y_continuous(
    limits = c(-1, 5),
    name   = "Treatment Effect"
  ) +
  ggtitle("Varying Distance to Turkish Coast") +
  theme_bw()
# Show the figure
figs710 Example AI Transparency Statement
Use of AI for Writing Support
Artificial intelligence tools were employed in a limited capacity to assist with the refinement of the written text. Specifically, I used AI to correct grammatical and syntactical errors in order to enhance readability and stylistic consistency. These interventions were strictly editorial in nature. The structure of the paper, the central argument, and the review of the relevant literature were developed independently in consultation with the advisor.
Use of AI for Coding Support
AI was also used to resolve syntactical errors in R code. This primarily concerned the specification of regression models and the formatting of graphical output. The role of AI in this context was limited to correcting technical implementation issues, ensuring that the code executed properly. All substantive decisions regarding model choice, robustness checks, coding strategy, and interpretation of results were determined independently.
Intellectual Contribution and Ownership
The intellectual contribution of this paper—including its theoretical framing, methodological design, and empirical analysis—rests entirely with the author. AI assistance was confined to technical and editorial functions, comparable to the support that might otherwise be obtained from proofreading services or coding documentation. I remain solely responsible for the scholarly content, interpretation of findings, and framing the contribution to the literature.