# =====================================================================
# Age-Based Screening for Driving Licence Renewal: Figures
#
# Replication code for: https://doi.org/10.1016/j.tranpol.2026.104493
#
# Produces four figures:
#   Fig. 1  Use of ABS in the European Union as of 2026, based on own research
#   Fig. 3  results.pdf            - Given gamma= 0: Disaggregated NPV of performing ABS for one year
#   Fig. 4  MCS.pdf                - Histogram of Monte Carlo Simulation results (n = 10,000)
#   Fig. 2  NPV_gamma.pdf          - Overall NPV of performing ABS for one year, depending on parameter choice for gamma
#
# Required input files (place in the working directory):
#   implementation.csv   semicolon-separated; columns: country, ISO-2, has_policy
#   results.csv          semicolon-separated; columns: effect, value
#   mcs.csv              semicolon-separated; column:  NPV
#   npv_gamma.csv        semicolon-separated; columns: gamma, npv
#
# Please set the working directory to the directory where the data files are stored.
#
# Fonts: To replicate the exact plot font, please install Latin Modern Roman 
# (lmroman10-regular.otf / lmroman10-bold.otf) on your system. Set "USE_CUSTOM_FONT"
# to TRUE and insert the paths in Section 1 to the location of these files.
#
# Tested with R 4.6.1, tidyverse 2.0.0
# =====================================================================

rm(list = ls())

library(tidyverse)
library(showtext)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
library(scales)
library(ggrepel)
library(ggh4x)


# =====================================================================
# 1. CONFIGURATION
# =====================================================================

# --- Paths -----------------------------------------------------------
path_policy <- "implementation.csv"
path_results <- "results.csv"
path_mcs     <- "mcs.csv"
path_npv     <- "npv_gamma.csv"
path_output  <- "."

# --- Font registration -----------------------------------------------
USE_CUSTOM_FONT <- F

font_regular <- ""
font_bold    <- ""

if (USE_CUSTOM_FONT && file.exists(font_regular) && file.exists(font_bold)) {
  font_add("LM Roman", regular = font_regular, bold = font_bold)
  plot_font <- "LM Roman"
} else {
  warning("Latin Modern font files not found; falling back to default device font.")
  plot_font <- ""
}

showtext_auto()

# --- Shared colours --------------------------------------------------
col_screening    <- "#ef7d00"   # orange: age-based screening in place
col_no_screening <- "#65b32e"   # green:  no age-based screening
col_non_eu       <- "grey90"
col_primary      <- "#00305D"   # dark blue for histogram and scatter


# =====================================================================
# 2. LOAD DATA
# =====================================================================

# --- Policy status by EU member state (own research, April 2026) -----
# has_policy: 1 = age-based screening required, 0 = not required
policy_data <- read.csv(path_policy, sep = ";") %>%
  mutate(
    has_policy = case_when(
      has_policy == 1 ~ "Age-based Screening",
      has_policy == 0 ~ "No Age-based Screening",
      TRUE            ~ NA_character_
    )
  )

res <- read.csv(path_results, sep = ";") %>%
  mutate(
    value  = as.numeric(gsub(",", "", value)),
    effect = factor(effect, levels = rev(unique(effect)))
  )

# --- Monte Carlo draws of the NPV ------------------------------------
mcs <- read.csv(path_mcs, sep = ";")

# --- NPV as a function of gamma --------------------------------------
npv <- read.csv(path_npv, sep = ";")


# =====================================================================
# 3. FIGURE 1 - POLICY STATUS ACROSS EU MEMBER STATES
# =====================================================================

# ---------------------------------------------------------------------
# 3.1 Join policy status onto country polygons
# ---------------------------------------------------------------------

world_sf <- ne_countries(scale = "medium", returnclass = "sf")

map_sf <- world_sf %>%
  mutate(
    has_policy = case_when(
      name %in% policy_data$country & name != "Russia" ~
        policy_data$has_policy[match(name, policy_data$country)],
      TRUE ~ "Non-EU"
    ),
    iso_label = case_when(
      name %in% policy_data$country ~
        policy_data$`ISO.2`[match(name, policy_data$country)],
      TRUE ~ NA_character_
    )
  )

# ---------------------------------------------------------------------
# 3.2 Label positions
# ---------------------------------------------------------------------

map_labels <- map_sf %>%
  filter(!is.na(iso_label) & iso_label != "Non-EU") %>%
  mutate(
    centroid = st_centroid(geometry, of_largest_polygon = TRUE),
    x = st_coordinates(centroid)[, 1],
    y = st_coordinates(centroid)[, 2]
  ) %>%
  st_drop_geometry()

map_labels_adj <- map_labels %>%
  rename(x_orig = x) %>%
  mutate(
    x = case_when(
      name == "Luxembourg" ~ x_orig + 1.4,
      name == "Malta"      ~ x_orig + 1.4,
      name == "Cyprus"     ~ x_orig - 1.8,
      name == "Greece"     ~ x_orig - 0.8,
      name == "Sweden"     ~ x_orig - 1.0,
      name == "Croatia"    ~ x_orig + 0.2,
      TRUE                 ~ x_orig
    ),
    # Only Croatia requires a vertical shift.
    y = case_when(
      name == "Croatia" ~ y + 0.4,
      TRUE              ~ y
    )
  )

# ---------------------------------------------------------------------
# 3.3 Leader lines
# ---------------------------------------------------------------------

bar_data <- map_labels_adj %>%
  filter(name %in% c("Luxembourg")) %>%
  mutate(
    x_anchor = x_orig,
    x_tip = case_when(
      name == "Luxembourg" ~ x - 0.9
    )
  )

# ---------------------------------------------------------------------
# 3.4 Plot
# ---------------------------------------------------------------------

p_map <- ggplot(map_sf) +
  geom_sf(aes(fill = has_policy), colour = "grey40", size = 0.15) +
  geom_segment(
    data = bar_data,
    aes(x = x_anchor, xend = x_tip, y = y, yend = y),
    colour = "black",
    linewidth = 0.3
  ) +
  geom_text(
    data = map_labels_adj,
    aes(x = x, y = y, label = iso_label),
    family = plot_font,
    size = 3,
    fontface = "bold"
  ) +
  scale_fill_manual(
    name = "Status",
    values = c(
      "Non-EU"                 = col_non_eu,
      "Age-based Screening"    = col_screening,
      "No Age-based Screening" = col_no_screening
    ),
    na.value = "lightgrey",
    guide = guide_legend(nrow = 3)
  ) +
  labs(
    title    = "ABS implementation in EU member states",
    subtitle = "As of April 2026",
    fill     = "Status"
  ) +
  coord_sf(xlim = c(-12, 35), ylim = c(33, 73), expand = FALSE) +
  theme_bw(base_size = 16) +
  theme(
    text            = element_text(family = plot_font),
    legend.position = "right",
    legend.title    = element_text(face = "bold"),
    axis.text       = element_blank(),
    axis.ticks      = element_blank(),
    axis.title      = element_blank(),
    panel.grid      = element_blank(),
    plot.title      = element_text(face = "bold", size = 18),
    plot.subtitle   = element_text(colour = "grey30"),
    plot.caption    = element_text(size = 8)
  )

ggsave(file.path(path_output, "Implementation_EU.pdf"),
       plot = p_map, width = 8, height = 5)


# =====================================================================
# 4. FIGURE 2 - COST AND BENEFIT COMPONENTS
# =====================================================================

# ---------------------------------------------------------------------
# 4.1 Component labels (plotmath expressions)
# ---------------------------------------------------------------------

label_map <- c(
  "GP employment"     = "SOC[me]",
  "Admin employment"  = "SOC[ad]",
  "Driver effort"     = "SOC[ti]",
  "Fee"               = "SOC[fe]",
  "Δ SOC"             = "bold('Δ')~bold(SOC)",
  
  "SS_cs"             = "SS[cs]",
  "SS_fc"             = "SS[fc]",
  "Δ SS"              = "bold('Δ')~bold(SS)",
  
  "EEs driving"       = "EE[d]",
  "EEs mental health" = "EE[de]",
  "EE Sm"             = "EE[su]",
  "Δ EEs"             = "bold('Δ')~bold(EE)",
  
  "Fuel tax"          = "GOV[fu]",
  "Parking VAT"       = "GOV[p]",
  "Car tax"           = "GOV[vt]",
  "GOV sm"            = "GOV[su]",
  "Δ GOV"             = "bold('Δ')~bold(GOV)",
  
  "NPV"               = "bold(NPV)"
)

res <- res %>%
  mutate(effect_char = as.character(effect)) %>%
  mutate(label_text = sapply(effect_char, function(x) {
    if (!is.na(label_map[x])) {
      label_map[x]
    } else {
      quoted <- gsub("'", "\\\\'", x)
      paste0("'", quoted, "'")
    }
  }))

# ---------------------------------------------------------------------
# 4.2 Fill gradient
# ---------------------------------------------------------------------

min_val  <- -1000e6
med_mid  <-  -200e6
low_mid  <-   -50e6
high_mid <-   100e6
max_val  <-   500e6

stops  <- c(min_val, med_mid, low_mid, high_mid, max_val)
values <- (stops - min_val) / (max_val - min_val)
colors <- c("red", "#F7766A", "#FFA48F", "#ADFFCA", "#00FF5A")

# ---------------------------------------------------------------------
# 4.3 Separator lines between the four effect blocks
# ---------------------------------------------------------------------

spine_ticks <- c("Δ SOC", "Δ SS", "Δ EEs", "Δ GOV")

ticks_df <- tibble(effect = spine_ticks) %>%
  mutate(
    effect = factor(effect, levels = levels(res$effect)),
    y      = as.numeric(effect) - 0.5
  )

# ---------------------------------------------------------------------
# 4.4 Value labels
# ---------------------------------------------------------------------

euro_mio <- label_number(
  scale        = 1e-6,
  suffix       = " Mio",
  prefix       = "€",
  big.mark     = ",",
  decimal.mark = ".",
  accuracy     = 0.1
)

res <- res %>%
  mutate(
    label_value = ifelse(
      value < 0,
      paste0("\u2212", euro_mio(abs(value))),
      euro_mio(value)
    )
  )

# ---------------------------------------------------------------------
# 4.5 Layout constants
# ---------------------------------------------------------------------

x_lim        <- max(abs(res$value), max_val) * 1.6
gap          <- x_lim * 0.11
x_range      <- c(-x_lim, x_lim)
label_offset <- x_lim * 0.01

# ---------------------------------------------------------------------
# 4.6 Plot
# ---------------------------------------------------------------------

p_results <- ggplot(res, aes(y = effect)) +
  # Bars left of the spine (negative values)
  geom_col(
    data = filter(res, value < 0),
    aes(x = value, fill = value, alpha = abs(value)),
    colour = "black", width = 0.6,
    position = position_nudge(x = -gap)
  ) +
  # Bars right of the spine (positive values)
  geom_col(
    data = filter(res, value >= 0),
    aes(x = value, fill = value, alpha = abs(value)),
    colour = "black", width = 0.6,
    position = position_nudge(x = gap)
  ) +
  # Block separators
  geom_hline(
    data = ticks_df,
    aes(yintercept = y),
    colour = "black", linewidth = 0.6
  ) +
  # Component labels in the central column
  geom_text(
    aes(x = 0, label = label_text),
    parse = TRUE,
    family = plot_font,
    size = 6
  ) +
  # Value labels, right-hand bars
  geom_text(
    data = filter(res, value >= 0),
    aes(x = value + gap + label_offset, label = label_value),
    hjust = 0, size = 6, family = plot_font
  ) +
  # Value labels, left-hand bars
  geom_text(
    data = filter(res, value < 0),
    aes(x = value - gap - label_offset, label = label_value),
    hjust = 1, size = 6, family = plot_font
  ) +
  scale_fill_gradientn(
    colours = colors,
    values  = values,
    limits  = c(min_val, max_val),
    oob     = scales::squish,
    guide   = "none"
  ) +
  # The two vertical spines bounding the label column
  geom_vline(
    xintercept = c(-gap, gap),
    colour = "black", linetype = "solid", linewidth = 0.5
  ) +
  scale_alpha_continuous(range = c(0.3, 1), guide = "none") +
  coord_cartesian(xlim = x_range) +
  scale_x_continuous(labels = NULL) +
  scale_y_discrete(limits = levels(res$effect), labels = NULL,
                   expand = c(0, 0.5)) +
  labs(x = NULL, y = NULL) +
  theme_bw(base_size = 14) +
  theme(
    text         = element_text(family = plot_font),
    panel.grid   = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_blank()
  )

ggsave(file.path(path_output, "results.pdf"),
       plot = p_results, width = 11.1, height = 10)


# =====================================================================
# 5. FIGURE 3 - MONTE CARLO DISTRIBUTION OF THE NPV
# =====================================================================

p_mcs <- ggplot(mcs, aes(x = NPV)) +
  geom_histogram(colour = "white", fill = col_primary, alpha = 0.7) +
  labs(x = "NPV in Million Euros", y = "Count") +
  scale_y_continuous(
    breaks = pretty_breaks(n = 3),
    expand = expansion(mult = c(0, 0.05))
  ) +
  scale_x_continuous(
    limits = c(-2100e6, 1500e6),
    labels = function(x) x / 1e6
  ) +
  theme_bw(base_size = 16) +
  theme(
    panel.grid.major.y = element_line(linewidth = 0.5, colour = "grey80"),
    panel.grid.minor   = element_blank(),
    axis.text          = element_text(size = 13, colour = "black"),
    axis.title         = element_text(size = 13, face = "bold"),
    text               = element_text(family = plot_font)
  )

ggsave(file.path(path_output, "MCS.pdf"),
       plot = p_mcs, width = 9, height = 3.5)


# =====================================================================
# 6. FIGURE 4 - NPV AS A FUNCTION OF GAMMA
# =====================================================================

p_npv <- ggplot(npv, aes(x = gamma, y = npv)) +
  geom_point(size = 0.01, colour = col_primary) +
  labs(x = NULL, y = NULL) +
  
  # Axis titles placed manually near the axis tips
  annotate(
    "text", x = -0.1, y = 1.5e9, label = "NPV [€]",
    family = plot_font, fontface = "bold", size = 4.6, hjust = 1
  ) +
  annotate(
    "text", x = 3.05, y = -5e8, label = "Gamma [%]",
    family = plot_font, fontface = "bold", size = 4.6, hjust = 1
  ) +
  
  scale_y_continuous(
    position     = "right",
    limits       = c(-2.5e9, 1.5e9),
    breaks       = seq(-5e9, 5e9, by = 5e8),   # major ticks every 500 m
    minor_breaks = seq(-5e9, 5e9, by = 1e8),   # minor ticks every 100 m
    guide        = guide_axis(minor.ticks = TRUE),
    labels = function(x) {
      # Label only multiples of 1e9, and never zero
      ifelse(x %% 1e9 == 0 & x != 0,
             trimws(format(x, scientific = TRUE)),
             "")
    },
    expand = expansion(mult = c(0, 0.05))
  ) +
  scale_x_continuous(
    limits       = c(-3.1, 3.1),
    expand       = c(0, 0),                    # no padding beyond the limits
    breaks       = seq(-4, 4, by = 0.5),       # major ticks every 0.5
    minor_breaks = seq(-4, 4, by = 0.1),       # minor ticks every 0.1
    guide        = guide_axis(minor.ticks = TRUE),
    labels = function(x) {
      # Label only integers; pad zero away from the vertical axis
      ifelse(x %% 1 == 0,
             ifelse(x == 0, "   0", as.character(x)),
             "")
    }
  ) +
  coord_axes_inside(labels_inside = TRUE) +
  theme_bw(base_size = 16) +
  theme(
    panel.grid   = element_blank(),
    panel.border = element_blank(),
    axis.line    = element_line(linewidth = 0.5, colour = "black"),
    axis.text    = element_text(size = 13, colour = "black"),
    text         = element_text(family = plot_font),
    
    # Ticks are thinner than the axis lines; minor ticks are half the
    # length of major ticks.
    axis.ticks              = element_line(linewidth = 0.25, colour = "black"),
    axis.ticks.length       = unit(0.2, "cm"),
    axis.minor.ticks.length = rel(0.5)
  )

ggsave(file.path(path_output, "NPV_gamma.pdf"),
       plot = p_npv, width = 8, height = 5)

# =====================================================================
# END OF SCRIPT
# =====================================================================