##################################################
################ PREPARATION #####################
##################################################
# Clear environment
rm(list=setdiff(ls(), "path_to"))

# Load data
world = read_dta(path_to("world_raw"))

# Load aux data
pop = read.csv(path_to("wdi_pop_total"), skip=4)
pop_share14 = read.csv(path_to("wdi_pop_share_14"), skip=4)
gdp = read.csv(path_to("wdi_gdp"), skip=4)
qogcs = read_dta(path_to("qogcs"))[, c("ccodealp", "wdi_lfpfilo15", "wdi_lfpmilo15", "gggi_ggi", "wdi_pop", "wdi_pop14", "wef_ler", "gii_gii"),]
romano2021 = read_excel(path_to("romano2021"), sheet = "Table 1", na="NA")
covid19 = read.csv(path_to("covid19"), as.is=T) %>% rename_with(tolower) %>% dplyr::select(
  iso3, total_vaccinations_per100, persons_vaccinated_1plus_dose_per100,
  persons_last_dose_per100, persons_booster_add_dose_per100
)
schulz2019 = read_dta(path_to("schulz2019")) %>% dplyr::select(js_code, BloodEstPopVol)
cohn2019 = read_dta(path_to("cohn2019"))


##################################################
################ CLEANING ########################
##################################################
# Merge QOG
qogcs = rename_with(qogcs, function(x) {paste0("qogcs_", x)})
world[!world$iso3 %in% qogcs$qogcs_ccodealp, c("iso3", "country")] # Missing countries
world = left_join(world, qogcs, by=c("iso3"="qogcs_ccodealp"))

# Merge GDP from WDI (pre-processed)
gdp = gdp %>% transmute(
    iso3 = Country.Code,
    gdp_pc_ppp_2010_2019 = rowMeans(pick(X2010:X2019), na.rm = TRUE),
    ln_gdp_pc_ppp_2010_2019 = log(gdp_pc_ppp_2010_2019)
)
gdp[gdp$iso3 == "VEN", "gdp_pc_ppp_2010_2019"] = 15937.53
gdp[gdp$iso3 == "VEN", "ln_gdp_pc_ppp_2010_2019"] = log(15937.53)
gdp = gdp %>% add_row(iso3="TWN", gdp_pc_ppp_2010_2019 = 45698.71, ln_gdp_pc_ppp_2010_2019 = log(45698.71))
world = world %>% left_join(gdp, by="iso3")
var_label(world$gdp_pc_ppp_2010_2019)    = "Average annual GDP 2010-2019, per capita, PPP-adjusted"
var_label(world$ln_gdp_pc_ppp_2010_2019) = "Log. average annual GDP 2010-2019, per capita, PPP-adjusted"

# Merge Romano et al. 2021 data
romano2021_means = romano2021 %>%
  group_by(Nations) %>%
  summarise(romano_coop = mean(`Mean_coop*10`, na.rm = TRUE))
romano2021_means$romano_coop_std = scale(romano2021_means$romano_coop) %>% as.numeric()
romano2021_means$iso3 = countrycode(romano2021_means$Nations, origin="country.name", destination="iso3c")
romano2021_means$Nations = NULL
world = world %>% left_join(romano2021_means, by="iso3")
var_label(world$romano_coop)     = "Cooperation measure from Romano et al. 2021b, country means, NOT standardized"
var_label(world$romano_coop_std) = "Cooperation measure from Romano et al. 2021b, std. country means"

# Merge Covid 19
covid19 = covid19 %>% rename(
    covid_total_vaccinations_pc = total_vaccinations_per100,
    covid_1plus_dose_pc         = persons_vaccinated_1plus_dose_per100,
    covid_fully_vaccinated_pc   = persons_last_dose_per100,
    covid_with_booster_pc       = persons_booster_add_dose_per100
  )
world = world %>% left_join(covid19, by = "iso3")
var_label(world$covid_total_vaccinations_pc) = "Total vaccinations (pc)"
var_label(world$covid_1plus_dose_pc)         = "At least one dose (%)"
var_label(world$covid_fully_vaccinated_pc)   = "Fully vaccinated (%)"
var_label(world$covid_with_booster_pc)       = "At least one booster (%)"

# Merge Schulz et al. (2019) data
world = world %>% left_join(schulz2019, by = join_by(iso3 == js_code))
world = world %>% dplyr::rename(c("schulz_BloodEstPopVol" = "BloodEstPopVol"))
world$schulz_BloodEstPopVol_std = scale(world$schulz_BloodEstPopVol) %>% as.numeric
var_label(world$schulz_BloodEstPopVol)     = "Blood donations per 1000 inhabitans (voluntary to non-family, WHO), data from Schulz et al. 2019"
var_label(world$schulz_BloodEstPopVol_std) = "Standardized blood donations per 1000 inhabitans (schulz_BloodEstPopVol)"

# Merge Cohn et al (2019) data
cohn2019 = cohn2019 %>% filter(cond %in% c(0, 1)) %>%
  mutate(cond = ifelse(cond == 0, "cohn_nomoney", "cohn_money")) %>%
  mutate(iso3 = countrycode(Country, origin = "country.name.en", destination = "iso3c")) %>%
  group_by(iso3, cond) %>%
  summarise(value = mean(response, na.rm = TRUE)) %>% ungroup %>%
  pivot_wider(names_from = cond, values_from = value) %>%
  mutate(cohn_rank = rank(cohn_nomoney, ties.method = "first"))
var_label(cohn2019$cohn_money) = "Cohn et al. (2019) - Money return rate"
var_label(cohn2019$cohn_nomoney) = "Cohn et al. (2019) - No money return rate"
var_label(cohn2019$cohn_rank) = "Cohn et al. (2019) - Rank in no money condition"
world = world %>% left_join(cohn2019, by = "iso3")

# Add population for missing countries
# Manually add WDI 2019 data that is missing for HKG and XKX in QOG
# XKX pop share below 15: Assume identical to SRB.
# TWN population: https://eng.stat.gov.tw/point.asp?index=9
# TWN pop share below 15: Assume identical to HKG, which seems to be an accurate
# assumption (https://www.statista.com/statistics/321439/taiwan-population-distribution-by-age-group/)
world[world$iso3 == "HKG", "qogcs_wdi_pop"] = pop[pop$Country.Code == "HKG", "X2019"]
world[world$iso3 == "XKX", "qogcs_wdi_pop"] = pop[pop$Country.Code == "XKX", "X2019"]
world[world$iso3 == "TWN", "qogcs_wdi_pop"] = 23600000
world[world$iso3 == "HKG", "qogcs_wdi_pop14"] = pop_share14[pop_share14$Country.Code == "HKG", "X2019"]
world[world$iso3 == "XKX", "qogcs_wdi_pop14"] = pop_share14[pop_share14$Country.Code == "SRB", "X2019"]
world[world$iso3 == "TWN", "qogcs_wdi_pop14"] = pop_share14[pop_share14$Country.Code == "HKG", "X2019"]

# Derive population 15+
world$pop15 = world$qogcs_wdi_pop * (1-world$qogcs_wdi_pop14/100)
var_label(world$pop15) = "Population (15+), total, based on WDI (qogcs) data"
if (any(is.na(world$pop15))) {message(
  paste("pop15 missing for:", paste0(world$iso3[is.na(world$pop15)], collapse=", "))
)}
world$qogcs_wdi_pop = world$qogcs_wdi_pop14 = NULL

# Derive weights, population 15
world$wgt_pop15 = world$pop15 / sum(world$pop15, na.rm=T)
var_label(world$wgt_pop15) = "Country weight, proportional to Population aged 15+, pop15/(sum of pop15)"


##################################################
################ EXPORT ##########################
##################################################
# Test for duplicate countries
if (any(duplicated(world$iso3))) {message(paste0("Duplicate countries: ", world$iso3[duplicated(world$iso3)]))}
world = world[!duplicated(world$iso3),]

# Save
saveRDS(world, path_to("world"))

