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/package" with the actual path to your package.
root_dir <- "your/path/to/package/"
root_dir <- "/Users/hectorsegura/Library/CloudStorage/OneDrive-Universidaddelosandes/World_Bank/AI_exposure/Reproducibility_package"
cleaned_data_dir <- file.path(root_dir, "Data", "Cleaned")
graphs_dir <- file.path(root_dir, "Outputs")
```{r Fig11, results = "show", fig.height=8, fig.width= 12}
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
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")
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/package" with the actual path to your package.
root_dir <- "your/path/to/package/"
root_dir <- "/Users/hectorsegura/Library/CloudStorage/OneDrive-Universidaddelosandes/World_Bank/AI_exposure/Reproducibility_package"
cleaned_data_dir <- file.path(root_dir, "Data", "Cleaned")
graphs_dir <- file.path(root_dir, "Outputs")
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
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")
GDP_pc_data <- WDI(country = "all", indicator = c("ny_gdp_pcap_pp_kd" = "NY.GDP.PCAP.PP.KD"),
latest = 1, extra = F)
GDP_pc_data <- WDI(country = "all", indicator = c("ny_gdp_pcap_pp_kd" = "NY.GDP.PCAP.PP.KD"),
latest = 1, extra = F)
GDP_pc_data <- GDP_pc_data %>%
rename(year_GDP = year, countrycode = iso3c) %>%
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))
