Disruption without Dividend?

How the Digital Divide and Task Differences split GenAI’s Global Impact

knitr::opts_chunk$set(comment = "#>",
                      echo = T,
                      results = 'show',
                      message = FALSE, 
                      warning = FALSE,
                      cache = F)

knitr::opts_knit$set(root.dir = here::here())

library(readxl)
library(openxlsx)
library(tidyverse)
library(wbstats)
library(WDI)
library(gghighlight)
library(lubridate)
library(ggridges)
library(sjPlot)
library(webshot)
library(DT)
library(tidytext)
library(plotly)
library(stringr)
library(stringdist)
library(fuzzyjoin)
library(wbstats)
library(countrycode)
library(dplyr)
library(tidyr)
library(purrr)
library(mgcv)
library(emmeans)
library(ggplot2)
library(flextable)
library(haven)
library(mgcv)
library(emmeans)
library(tibble)

my_colors <- c("Gradient 1" = "steelblue", "Gradient 2" = "skyblue", 
               "Gradient 3" = "orange2", "Gradient 4" = "#D7191C", 
               "Minimal Exposure" = "lightgreen",  "Not Exposed" = "#BDBDBD")
my_shapes <- c("Gradient 1" = 15, "Gradient 2" = 16, 
               "Gradient 3" = 18, "Gradient 4" = 17, 
               "Minimal Exposure" = 19, "Not Exposed" = 19) 
my_sizes <- c("Gradient 1" = 3, "Gradient 2" = 3, 
              "Gradient 3" = 4, "Gradient 4" = 4, 
              "Minimal Exposure" = 2, "Not Exposed" = 2) 
my_alpha <- c("Gradient 1" = 0.8, "Gradient 2" = 0.8, 
              "Gradient 3" = 0.8, "Gradient 4" = 0.8, 
              "Minimal Exposure" = 0.8, "Not Exposed" = 0.8) 

# Replace "your/path/to/Reproducibility_package" with the actual path to your package.

root_dir <- "C:/Users/wb519501/OneDrive - WBG/Desktop/RR_WLD_2025_528/files upload/v3/Reproducibility_package"

cleaned_data_dir <- file.path(root_dir, "Data", "Cleaned")

graphs_dir <- file.path(root_dir, "Outputs")

Figure 9

welfare <- read_excel(file.path(cleaned_data_dir, "Welfare", "Welfare_at_2d.xlsx")) %>% 
  group_by(countrycode, countryname) %>% 
  mutate(median_welfare = median(welfare),
         mean_welfare = mean(welfare)) %>% 
  ungroup() %>% 
  mutate(distance_median_welfare = (welfare - median_welfare)/median_welfare) %>% 
  mutate(ISCO_2d = as.double(ISCO_2d)) 

labels2d <- read_csv(file.path(dirname(cleaned_data_dir), "Raw", "Other", "mapping.csv")) %>% 
  select(isco08_2d, description_2d) %>% 
  distinct() %>% 
  mutate(label2d = paste(isco08_2d, "-", description_2d)) %>% 
  rename(ISCO_2d = isco08_2d) %>% 
  select(1,3)%>% 
  mutate(ISCO_2d = as.double(ISCO_2d)) 

welfare2 <- left_join(welfare, labels2d) %>% 
  mutate(ISCO_2d = factor(ISCO_2d),
         labels2d = factor(label2d)) %>%
  group_by(countryname, countrycode) %>% 
  mutate(total_empl = sum(weight)) %>% 
  ungroup() %>% 
  pivot_longer(cols = Gradient_1:Not_exposed,
               names_to = "potential25",
               values_to = "exposure_value") %>% 
  mutate(weight_by_exposure = weight * exposure_value) %>% 
  mutate(exposure_share_empl = weight_by_exposure / total_empl * 100) %>% 
  group_by(countryname, countrycode) %>% 
  mutate(check = sum(exposure_share_empl)) %>% 
  ungroup() %>% 
  mutate(welfare_distance_filtered = ifelse(exposure_value > 0 &
                                              !potential25 %in% c("Low_exposure", 
                                                                 "Not_exposed"), distance_median_welfare, NA)) 

# Remove outliers
outlier <- welfare2 %>% filter(distance_median_welfare > 6)

# Define custom legend colors
custom_Legend_colors <- setNames(
  c("red", "orange", "skyblue2", "blue"),
  c("Gradient_4", "Gradient_3", "Gradient_2", "Gradient_1"))

set.seed(1234)

# Plot
ggplot(welfare2, aes(x = fct_rev(label2d), y = distance_median_welfare)) +
  geom_jitter(aes(size = exposure_share_empl), alpha = 0.3, color = "grey") + # Plot all 
  geom_jitter(aes(x = fct_rev(label2d), 
                  y = welfare_distance_filtered, 
                  color = potential25, 
                  size = exposure_share_empl),
              alpha = 0.8) + # Ovelay exposed jobs
  geom_smooth(method = "loess", se = FALSE, 
              color = "black", linewidth = 0.3, 
              linetype ="dashed", aes(group = potential25)) +
  geom_hline(yintercept = 0, color = "brown", size = 0.7) +
  scale_color_manual(values = custom_Legend_colors) +  # Apply custom color scale
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
        legend.position = "top",  
        legend.box = "horizontal", 
        legend.box.just = "left",  
        legend.direction = "horizontal",
        legend.title = element_text(size = 10),
        axis.title.y = element_text(size = 10)) +  
  labs(x = "",
       y = "Welfare: distance from the median welfare (as % of the median)",
       size = "Size: share of total employment",
       alpha = "Color and shading: Exposure Gradient") +
  guides(
    color = guide_legend(order = 3, title.position = "top", nrow = 1, byrow = TRUE, title = NULL),
    size = guide_legend(order = 1, title.position = "top", nrow = 1, byrow = TRUE),
    alpha = "none") +
  ylim(-1, 4)

# Save plot
ggsave(filename = file.path(graphs_dir, "Main/Figures/genAI_exposure_task_content", "Figure9.jpeg"),
       width = 8, height = 9,
       units = "in",
       dpi = 450)

# List countries in the plot
welfare_countries <- welfare2 %>% select(countryname) %>% distinct()

Figure 11

data4d <- read_dta(file.path(cleaned_data_dir, "PIAAC/PIAAC_4digit_with_4digit_AI_tasks_scores.dta")) %>%
  zap_label() %>% 
  pivot_longer(NRA_Belgium: last_col(),
               names_to = "var",
               values_to = "value") %>% 
  drop_na(value) %>% 
  separate(var, into = c("var", "country"), sep = "_", extra = "merge") %>% 
  mutate(total_obs = ifelse(var == "total", value, NA)) %>% 
  fill(total_obs, .direction = "up") %>% 
  filter(!var %in% c("total", "NRM", "RC")) %>% 
  mutate(concept = case_when(var == "CU" ~ "Computer Use",
                             var == "NRA" ~ "Non-Routine Analytical",
                             var == "NRI" ~ "Non-Routine Interpersonal",
                             var == "RM" ~ "Routine Manual")) %>% 
  mutate(label4d = paste(isco4c, "-", OccupationName)) %>% 
  select(-c(isco4c, OccupationName)) %>% 
  select(label4d, everything()) %>% 
  mutate(concept = factor(concept, levels = c("Computer Use",
                                              "Routine Manual",
                                              "Non-Routine Analytical",
                                              "Non-Routine Interpersonal"))) %>% 
  mutate(Exposure = ifelse(Exposure == "Low Exposure", "Minimal Exposure", Exposure),
         Exposure = str_trim(Exposure, side = "both")) %>% 
  mutate(label_outliers = paste(country, ":", label4d))

### Filter obs > 5
data_5obs <- data4d %>% filter(total_obs > 5) 

# List countries
countries_fig11 <- data_5obs %>% distinct(country)

### Mark outliers
outliers <- data_5obs %>% 
  filter( (concept == "Routine Manual" & value > -0.7 & Exposure == "Gradient 4") |
          (concept == "Routine Manual" & value > 0 & Exposure == "Gradient 3") |
          (concept == "Non-Routine Interpersonal" & value > 1 & Exposure == "Gradient 4") |
          (concept == "Routine Manual" & value > -0.7 & Exposure == "Gradient 4") |
          (concept == "Computer Use" & value < -0.3 & Exposure %in% c("Gradient 3", "Gradient 4"))
          )

ggplot(data_5obs, 
       aes(x = MeanScore2025, y = value)) + 
  geom_point(aes(color = Exposure,
                 shape = Exposure,
                 size = Exposure)) + 
  geom_smooth(method = "loess", se = FALSE) +  
  facet_wrap(~concept, scales = "free") + 
  theme(legend.position = "bottom") +
  labs(title = "",
       subtitle = "Countries with 4-digit PIAAC 2025 data: Belgium, Chechia, Chile, France, Poland, Slovak Rep., Spain",
       x = "AI exposure score, mean for occupation, ILO 2025)",
       y = "PIAAC task intensity value") +
  scale_color_manual(values = my_colors, 
                     breaks = c("Not Exposed", "Minimal Exposure", "Gradient 1", "Gradient 2", 
                                "Gradient 3", "Gradient 4")) +
  scale_shape_manual(values = my_shapes, 
                     breaks = c("Not Exposed", "Minimal Exposure", "Gradient 1", "Gradient 2", 
                                "Gradient 3", "Gradient 4")) +
  scale_size_manual(values = my_sizes, 
                    breaks = c("Not Exposed", "Minimal Exposure", "Gradient 1", "Gradient 2", 
                               "Gradient 3", "Gradient 4")) +
  theme(strip.text = element_text(size = 12),
        plot.title = element_text(size = 12),
        legend.title = element_blank(),
        legend.text = element_text(size = 12),
        axis.title = element_text(size = 12)
        ) +
  ggrepel::geom_label_repel(data = outliers, 
                            aes(x = MeanScore2025, y = value, label = label_outliers), 
                            size = 3, 
                            max.overlaps = 30,  # Allow more labels
                            box.padding = 0.7,  # Increase space around labels
                            segment.curvature = -0.2,  # Curve the connecting lines
                            segment.angle = 20)  # Angled connectors

# Save 
ggsave(filename = file.path(graphs_dir, "Main/Figures/genAI_exposure_task_content", "Figure11.png"),
       dpi = 400,
       height = 8,
       width = 12)

Figure 12

# Prepare data
data_clean <- data_5obs %>%
  mutate(concept = as.factor(concept),
         country = as.factor(country))

### Regressions: flexible curve for each country
models_country_smooth <- data_clean %>%
  group_split(concept) %>%
  set_names(map_chr(., ~ as.character(unique(.x$concept)))) %>%
  map(~ gam(value ~ s(MeanScore2025, country, bs = "fs"), data = .x))

data_country_smooth <- map2_dfr(
  models_country_smooth,
  names(models_country_smooth),
  ~ {.data <- .x$model
     .data$fit <- predict(.x, type = "response")
     .data$concept <- .y
     .data})

# Set exposure levels for margins
exposure_levels <- c(0.1, 0.25, 0.5, 0.75)  

# Create a prediction table at specified exposure levels
predicted_at_levels <- map2_dfr(
  models_country_smooth,
  names(models_country_smooth),
  ~ {
    model <- .x
    concept_name <- .y
    newdata <- expand.grid(
      MeanScore2025 = exposure_levels,
      country = levels(model$model$country)
    )
    preds <- predict(model, newdata = newdata, se.fit = TRUE)
    
    tibble(
      concept = concept_name,
      country = newdata$country,
      exposure = newdata$MeanScore2025,
      fit = preds$fit,
      se = preds$se.fit,
      ci_low = preds$fit - 1.96 * preds$se.fit,
      ci_high = preds$fit + 1.96 * preds$se.fit
    ) %>%
      mutate(across(c(fit, se, ci_low, ci_high), ~ round(.x, 3)))
  }
) %>% 
  mutate(concept = factor(concept, levels = c("Computer Use",
                                              "Routine Manual",
                                              "Non-Routine Analytical",
                                              "Non-Routine Interpersonal")))


ggplot(predicted_at_levels, aes(x = exposure, y = fit, color = country, group = country)) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.02) +
  facet_wrap(~ concept, scales = "free_y") +
  labs(
    title = "Predicted Task Intensity by AI Exposure Score and Country",
    subtitle = "At selected exposure levels (0.1, 0.25, 0.5, 0.75) from country-specific smooths",
    x = "GenAI Exposure Score (2025)",
    y = "Predicted Task Intensity"
  ) +
  theme(legend.position = "bottom")

# Save
ggsave(filename = file.path(graphs_dir, "Main/Figures/genAI_exposure_task_content", "Figure12.png"),
       dpi = 400,
       height = 8,
       width = 12)

Figure 13

GDP_pc_data <- readRDS(file.path(dirname(cleaned_data_dir), "Raw", "Other", "WDI_GDP_region_data_latest24.rds"))

GDP_pc_data <- GDP_pc_data %>% 
  rename(year_GDP = year, countrycode = iso3c, ny_gdp_pcap_pp_kd = NY.GDP.PCAP.PP.KD) %>% 
  select(countrycode, year_GDP, ny_gdp_pcap_pp_kd) %>% 
  zap_label()

country_region_data <- wb_countries()

country_region_data <- country_region_data %>% 
  rename(countrycode = iso3c, 
         region = region_iso3c,
         regionname = region, 
         adminregion = admin_region_iso3c,
         adminregionname = admin_region, 
         incomelevel = income_level_iso3c,
         incomelevelname = income_level) %>% 
  select(countrycode, region, regionname, adminregion, adminregionname, incomelevel, incomelevelname)

new_data_2d <- read_dta(file.path(cleaned_data_dir, "PIAAC/2d_step_piaac_for_pawel.dta")) %>% 
  zap_label()

new_data_2d <- new_data_2d %>%
  select(-starts_with("weight_")) %>% 
  pivot_longer(cols = c(starts_with("NRA_"), starts_with("NRI_"), starts_with("RC_"), 
                        starts_with("RM_"), starts_with("NRM_"), starts_with("CU_"),
                        starts_with("total_")),
               names_to = c(".value", "countryname"),
               names_pattern = "([^_]+_)(.*)") %>% 
  select(isco4c, Exposure, ISCO081digit, OccupationName, MeanScore2025, SD, isco2c, 
         NRA_, NRI_, RC_, RM_, NRM_, CU_, total_, countryname) %>% 
  mutate(countryname = countrycode(countryname, origin = "country.name", destination = "country.name",
                                   nomatch = countryname),
         countrycode = countrycode(countryname, origin = "country.name", destination = "iso3c", 
                                   custom_match = c("Kosovo" = "XKX"))) %>% 
  mutate(step = ifelse(countrycode %in% c("ARM", "BOL", "COL", "SLV", "GEO", "GHA", "KEN", "XKX", "LAO",
                                          "MKD", "PHL", "SRB", "LKA", "UKR", "VNM"), 1, 0), 
         .after = countryname)

new_data_2d <- new_data_2d %>% 
  left_join(country_region_data, by = "countrycode") %>% 
  left_join(GDP_pc_data %>% select(-year_GDP), by = "countrycode") %>% 
  arrange(countryname) %>% 
  filter(countryname != "Yunnan") %>% 
  mutate(Exposure = str_trim(Exposure, side = "both"))

countries_new_file <- new_data_2d %>% distinct(countrycode, countryname) # 46 countries 
codes_new_file <- new_data_2d %>% distinct(countrycode) # 46 countries 
names_new_file <- new_data_2d %>% distinct(countryname) # 46 countries 

# Process data
new_data_4d_clean <- new_data_2d %>% 
  filter(!is.na(isco4c)) %>% 
  mutate(isco4c = as.character(isco4c)) %>% 
  pivot_longer(NRA_:CU_,
               names_to = "concept",
               values_to = "value") %>% 
  drop_na(value) %>% 
  rename(total_obs = total_) %>% 
  arrange(countrycode, isco2c) %>% 
  rename(MeanScore4d = MeanScore2025,
         label1d = ISCO081digit) %>% 
  group_by(countrycode, isco2c) %>% 
  mutate(MeanScore2d = round(mean(MeanScore4d), 2)) %>% 
  ungroup() %>% 
  filter(!concept %in% c("total", "NRM_", "RC_")) %>% 
  mutate(concept = case_when(concept == "CU_" ~ "Computer Use",
                             concept == "NRA_" ~ "Non-Routine Analytical",
                             concept == "NRI_" ~ "Non-Routine Interpersonal",
                             concept == "RM_" ~ "Routine Manual")) %>% 
  select(isco2c, everything()) %>% 
  mutate(concept = factor(concept, levels = c("Computer Use",
                                              "Routine Manual",
                                              "Non-Routine Analytical",
                                              "Non-Routine Interpersonal"))) %>% 
  select(1:8, value, everything())  %>% 
  mutate(Exposure = ifelse(Exposure == "Low Exposure", "Minimal Exposure", Exposure)) 

# Create source mapping
country_names <- new_data_4d_clean %>% distinct(countryname) # 46 country names

# Assign source labels to country data
sources_df <- tibble(
  countryname = country_names$countryname,
  source = case_when(
    countryname %in% c("Bolivia", "Colombia", "Laos", "Sri Lanka", "Vietnam") ~ "STEP 2012",
    countryname %in% c("Armenia", "El Salvador", "Georgia", "Ghana", "Kenya", "North Macedonia", "Ukraine") ~ "STEP 2013",
    countryname %in% c("Serbia", "Kosovo", "Philippines") ~ "STEP 2015/16",
    countryname %in% c("Ecuador", "Greece", "Kazakhstan", "Mexico", "Peru", "Russia", "Slovenia", "Turkey") ~ "PIAAC first wave",
    countryname %in% c("Belgium", "Chile", "Czechia", "Denmark", "France", "Germany", "Hungary", "Ireland", "Israel",
                   "Italy", "Japan", "South Korea", "Lithuania", "Netherlands", "New Zealand", "Norway", "Poland",
                   "Singapore", "Slovakia", "Spain", "Sweden", "United Kingdom", "United States") ~ "PIAAC second wave",
    TRUE ~ NA_character_))

# Add sources to the main file - generate additional var STEP/PIAAC
step4d_full <- left_join(new_data_4d_clean, sources_df) %>% 
  mutate(source2 = ifelse(source %in% c("STEP 2012", "STEP 2013", "STEP 2015/16"), "STEP", "PIAAC"))

# LOAD separate file with employment weights at 2d
weights_data_2d <- read_dta(file.path(cleaned_data_dir, "PIAAC/tasks_exposure_2d_for_pawel.dta")) %>%
  zap_label() %>% 
  select(isco2c, countryname, weight) %>% 
  drop_na() %>% 
  rename(weight2d = weight) %>% 
  mutate(countryname = countrycode(countryname, origin = "country.name", destination = "country.name", 
                                   nomatch = countryname)) 

# Attach weights to the main dataset
new_data_weights_4d <- left_join(step4d_full, weights_data_2d)


# GAM WITH WEIGHTS for Gradients and Intensity
# 0) Ensure factors + weights
new_data_weights_4d_f <- new_data_weights_4d %>%
  mutate(incomelevel = factor(incomelevel),
    # Set ordering for Exposure
    Exposure = factor(Exposure, levels = c("Not Exposed", "Minimal Exposure", 
                                           "Gradient 1", "Gradient 2", 
                                           "Gradient 3", "Gradient 4"), 
                                            ordered = TRUE),
    weight = as.numeric(weight2d),
    weight2 = as.numeric(weight * ny_gdp_pcap_pp_kd)) %>%
    droplevels()

# 1) Fit one weighted model per concept (no smooth; Exposure is a factor)
models_exposure <- new_data_weights_4d_f %>%
  group_split(concept) %>%
  set_names(map_chr(., ~ as.character(unique(.x$concept)))) %>%
  map(~ gam(
    value ~ Exposure * incomelevel,  # categorical interaction
    data = .x,
    weights = .x$weight2))

# 2) Predictions at each Exposure x incomelevel combination (with SE & CI)
predicted_exposure <- map2_dfr(
  models_exposure,
  names(models_exposure),
  ~ {
    model <- .x
    concept_name <- .y
    newdata <- expand.grid(
      Exposure    = levels(model$model$Exposure),
      incomelevel = levels(model$model$incomelevel),
      KEEP.OUT.ATTRS = FALSE
    )
    preds <- predict(model, newdata = newdata, se.fit = TRUE)
    tibble(
      concept     = concept_name,
      Exposure    = newdata$Exposure,
      incomelevel = newdata$incomelevel,
      fit         = preds$fit,
      se          = preds$se.fit,
      ci_low      = preds$fit - 1.96 * preds$se.fit,
      ci_high     = preds$fit + 1.96 * preds$se.fit
      )})

# 3) Plot: lines across ordered Exposure categories, by income level, faceted by concept
ggplot(predicted_exposure %>%
    mutate(Exposure_num = as.numeric(Exposure), # to connect in order
           incomelevel = factor(incomelevel, levels = c("HIC", "UMC", "LMC")),
           concept = factor(concept, levels = c("Computer Use", "Routine Manual", 
                                               "Non-Routine Analytical", "Non-Routine Interpersonal"))),  
  aes(x = Exposure_num, y = fit, 
      color = incomelevel, group = incomelevel)) +
  geom_line(size = 0.7) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.1) +
  scale_x_continuous(
    breaks = seq_along(levels(predicted_exposure$Exposure)),
    labels = levels(predicted_exposure$Exposure)) +
  facet_wrap(~ concept, scales = "free_y") +
  labs(title = "",
       subtitle = "Weighted GAMs (parametric) with Exposure × incomelevel",
       x = "",
       y = "Predicted Task Intensity (PIAAC/STEP)") +
  theme(legend.position = "right", 
        legend.title = element_blank(),
        axis.text.x = element_text(angle = 90))

# Save plot
ggsave(filename = file.path(graphs_dir, "Main/Figures/genAI_exposure_task_content", "Figure13.png"),
       dpi = 400,
       height = 5,
       width = 8)

Figure 14

# Dataset: from elements in Figure 13
new_data_2d_clean1 <- new_data_4d_clean %>% 
  distinct(countrycode, countryname, label1d, isco2c, MeanScore2d, total_obs, value, concept, region, 
           regionname, adminregion, adminregionname, incomelevel, incomelevelname, ny_gdp_pcap_pp_kd) 

step2d_full <- left_join(new_data_2d_clean1, sources_df) %>% 
  mutate(source2 = ifelse(source %in% c("STEP 2012", "STEP 2013", "STEP 2015/16"), "STEP", "PIAAC"))

# LOAD separate file with employment weights at 2d
weights_data_2d <- read_dta(file.path(cleaned_data_dir, "PIAAC/tasks_exposure_2d_for_pawel.dta")) %>%
  zap_label() %>% 
  select(isco2c, countryname, weight) %>% 
  drop_na() %>% 
  rename(weight2d = weight) %>% 
  mutate(countryname = countrycode(countryname, origin = "country.name", destination = "country.name", 
                                   nomatch = countryname)) 
  

new_data_weights_2d <- left_join(step2d_full, weights_data_2d)

new_data_weights_2d_f <- new_data_weights_2d %>% 
  mutate(incomelevel = as.factor(incomelevel)) %>% 
  mutate(concept = recode(concept,
    "Computer Use" = "CU",
    "Routine Manual" = "RM",
    "Non-Routine Analytical" = "NRA",
    "Non-Routine Interpersonal" = "NRI"
  )) %>% 
  mutate(weight_alt = as.numeric(weight2d * ny_gdp_pcap_pp_kd))

# Create income and region labels
labels_country <- new_data_weights_2d_f %>% 
  select(countrycode, countryname, isco2c, label1d, region:source2) %>% 
  distinct() %>% 
  select(-c(isco2c, label1d)) %>% 
  distinct()

wide_data <- new_data_weights_2d_f %>%
  select(countrycode, countryname, incomelevel, isco2c, MeanScore2d, concept, value, weight2d, weight_alt) %>%
  pivot_wider(names_from = concept, values_from = value)

# GAM with k = 4
model_gam <- gam(
  MeanScore2d ~ s(CU, k = 4) + s(RM, k = 4) + s(NRA, k = 4) + s(NRI, k = 4), 
  data = wide_data,
  weights = weight2d)  # employment shares normalized within country

# Predict based on the model
wide_data$AdjustedExposure <- predict(model_gam, newdata = wide_data)

# Reshape to long format for plotting
data_country <- wide_data %>% 
  select(countryname, isco2c, MeanScore2d, AdjustedExposure, weight2d) %>% 
  group_by(countryname) %>% 
  summarise(mean_w = weighted.mean(MeanScore2d, weight2d, na.rm = TRUE), 
            adjusted_w = weighted.mean(AdjustedExposure, weight2d, na.rm = TRUE)) %>% 
  ungroup() %>% 
  left_join(labels_country) 

country_long <- data_country %>% 
  pivot_longer(cols = mean_w:adjusted_w,
               names_to = "type",
               values_to = "value") %>%
  mutate(type = recode(type,
                  mean_w = "Original Exposure",
                  adjusted_w = "Adjusted Exposure")) %>% 
  mutate(type = factor(type, levels = c("Original Exposure", "Adjusted Exposure"))) %>% 
  mutate(incomelevel = factor(incomelevel, levels = c("LMC", "UMC", "HIC")))

# Plot
ggplot(country_long,
       aes(x = log(ny_gdp_pcap_pp_kd), y = value)) +
  geom_point(aes(color = incomelevel)) +
  facet_wrap(~type) +
  geom_smooth(method = "gam", se = FALSE, color = "black") +  
  ggrepel::geom_text_repel(aes(label = countrycode)) +
  scale_color_manual(
    values = c(
      "HIC" = "red2",
      "UMC" = "blue",
      "LMC" = "green2")) +
  labs(x = "GDP per capita (log scale)",
       y = "Mean country-level exposure (weighted by 2-digit employment)") +
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        strip.text = element_text(size = 12),
        axis.title = element_text(size = 14))

# Save
ggsave(filename = file.path(graphs_dir, "Main/Figures/genAI_exposure_task_content", "Figure14.png"),
       dpi = 400,
       height = 8,
       width = 10)

Figure A15

data <- readRDS(file.path(dirname(cleaned_data_dir), "Raw", "Other", "GDP_internet_from_WDI.rds"))

data1 <- data %>%  
  # Drop rows with missing values in the selected indicators
  filter(!is.na(NY.GDP.PCAP.PP.KD) & !is.na(IT.NET.USER.ZS)) %>% 
  # Drop aggregates and regional NAs
  filter(!is.na(region),
         region != "Aggregates") %>% 
  # Select the latest year for each country
  group_by(country) %>%
  filter(year == max(year)) %>%
  ungroup() %>% 
  # Drop rows with missing year
  filter(!is.na(year)) %>% 
  # Filter data for the years 2019 or later
  filter(year >= 2019) %>% 
  # Create new variable lgdp (log of GDP per capita) and LAC
  mutate(lgdp = log(NY.GDP.PCAP.PP.KD)) %>% 
  # Create LAC identifier
  mutate(lac = ifelse(region =="Latin America & Caribbean", "LAC", "non-LAC"))

# Plot
ggplot(data1, aes(x = lgdp, y = IT.NET.USER.ZS, color = lac)) +
  geom_point(size = 1) +
  geom_smooth(method = "loess", formula = y ~ log(x), se = FALSE, aes(color = "Fitted line (loess)"), linewidth = 1) +
  scale_color_manual(values = c("Fitted line (loess)" = "green", "LAC" = "red", "non-LAC" = "blue"),
                     labels = c("Fitted line (loess)", "LAC", "non-LAC"),
                     guide = guide_legend(reverse = TRUE)) +
  labs(x = "GDP per capita, USD 2021 PPP (log scale) (WDI, latest year)", 
       y = "Internet Users per 100 people (ITU, latest, mostly 2023/24)", 
       color = "") +
  theme_minimal() +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 10),
        axis.title = element_text(size = 11))

# Save
ggsave(filename = file.path(graphs_dir, "Annex/Figures", "FigureA15_internet.jpeg"),
       width = 8, height = 6,
       units = "in",
       dpi = 300)