OHI Science | Citation policy +
+
OHI Science | Citation policy
+
This script downloads WGI data and prepares it for a pressures (1 - +WGI) and resilience data layer.
+No methods updates; additional year of data added.
+Potential improvement for future assessments: Explore a better way to +gapfill… based on a linear model (~line 166)
+Reference: http://info.worldbank.org/governance/wgi/index.aspx#home
+Downloaded:
+October 3, 2022 (source data updated Sep 2022)
+Description:
+The Worldwide Governance Indicators (WGI) project reports aggregate and
+individual governance indicators for 215 economies over the period
+1996–2021, for six dimensions of governance:
Time range: 1996-2021
+library(ohicore) # devtools::install_github('ohi-science/ohicore@dev')
+library(tools)
+library(dplyr)
+library(tidyr)
+library(WDI) # install.packages('WDI') # used to extract World Development Indicator (World Bank) data
+library(stringr)
+library(here)
+library(ggplot2)
+library(plotly)
+library(readr)
+
+
+# check website to see what years are available: http://info.worldbank.org/governance/wgi/index.aspx#home
+yr_start = 1996
+yr_end = 2021
+
+scen_year <- "2022"Download each of the 6 WGI indicators:
+## access data through the WDI package
+
+## get description of variables:
+indicators <- data.frame(WDI_data[[1]])
+indicators[grep("VA.EST", indicators$indicator), ]
+indicators[grep("PV.EST", indicators$indicator), ]
+indicators[grep("GE.EST", indicators$indicator), ]
+indicators[grep("RQ.EST", indicators$indicator), ]
+indicators[grep("RL.EST", indicators$indicator), ]
+indicators[grep("CC.EST", indicators$indicator), ]
+
+# identify the six indicators
+# WDIsearch('violence')# general search
+key_voice = WDI(
+ indicator = WDIsearch('Voice and Accountability: Estimate', field='name')['indicator'],
+ country = 'all', start = yr_start, end=yr_end)
+
+key_polst = WDI(
+ WDIsearch('Political Stability and Absence of Violence/Terrorism: Estimate', field='name')['indicator'],
+ country='all',start = yr_start, end=yr_end)
+
+key_gvtef = WDI(
+ WDIsearch('Government Effectiveness: Estimate', field='name')['indicator'],
+ country='all',start = yr_start, end=yr_end)
+
+key_regqt = WDI(
+ WDIsearch('Regulatory Quality: Estimate', field='name')['indicator'],
+ country='all',start = yr_start, end=yr_end)
+
+key_rolaw = WDI(
+ WDIsearch('Rule of Law: Estimate', field='name')['indicator'],
+ country='all',start = yr_start, end=yr_end)
+
+key_corrp = WDI(
+ WDIsearch('Control of Corruption: Estimate', field='name')['indicator'],
+ country='all',start = yr_start, end=yr_end)Combine the indicators into a single table, with a column for each +indicator, and rows for each country-year pair.
+d = key_voice %>%
+ select(country, year, indicator) %>%
+ left_join(key_polst %>% select(-iso2c), by=(c('country', 'year'))) %>%
+ left_join(key_gvtef %>% select(-iso2c), by=(c('country', 'year'))) %>%
+ left_join(key_regqt %>% select(-iso2c), by=(c('country', 'year'))) %>%
+ left_join(key_rolaw %>% select(-iso2c), by=(c('country', 'year'))) %>%
+ left_join(key_corrp %>% select(-iso2c), by=(c('country', 'year'))); head(d); summary(d); sapply(d, class)
+
+colnames(d) <- c("country", "year", "VA.EST", "PV.EST", "GE.EST", "RQ.EST", "RL.EST", "CC.EST")Uncomment the code chunk lines when updating WGI data, this will most +likely occur when calculating for new assessment year:
+date <- Sys.Date()
+write.csv(d, here(sprintf('globalprep/prs_res_wgi/v%s/raw/worldbank_wgi_from_wdi_api_%s.csv', scen_year, date)), row.names=FALSE)
+
+
+# This dataset currently has non-OHI regions includedThe first gapfilling occurs when we use the average of previous years +data within each region/indicator. This occurs when a region has data +for an indicator, but not for all years.
+Read in WGI data - change appended date in file name to +reflect the most recent version of the saved WGI data:
+d <- read.csv(here(paste0('globalprep/prs_res_wgi/v', scen_year, '/raw/worldbank_wgi_from_wdi_api_2022-10-03.csv'))) # change appended date to most recent versiond <- d %>%
+ pivot_longer(cols = VA.EST:CC.EST, names_to = "indicator", values_to = "value")
+
+
+## each country has 20 years of data
+d_gap_fill <- d %>%
+ group_by(country, year) %>%
+ mutate(NA_count_c_y = sum(is.na(value))) %>% # gf record: NA values within a region/year prior to gapfilling, max value is 6 (meaning that a country has no data)
+ ungroup() %>%
+ group_by(country, indicator) %>% # gapfill missing data with mean of values across years within the same region/indicator
+ mutate(ind_mean_c_i = mean(value, na.rm=TRUE)) %>%
+ ungroup() %>%
+ mutate(value = ifelse(is.na(value), ind_mean_c_i, value)) %>%
+ group_by(country, year) %>%
+ mutate(NA_count_post_gf1 = sum(is.na(value))) # gf record: NA values within a region/year after within region/indicator gapfilling (i.e. indicator is gapfilled by other years of data), used to cut regions <4 indicators (below) Once gapfilling is complete, the WGI scores are calculated as an +average of the 6 indicators. However, if a country is missing 4 or more +of the indicators within a year the average would be very biased. In +these cases, a different method should be used to gapfill these data
+(NOTE: for the 2022 assessment all regions had at least 4 of the 6 +indicators).
+countries_no_data <- d_gap_fill %>%
+ filter(NA_count_post_gf1 > 3)
+
+countries_no_data <- unique(countries_no_data$country)
+countries_no_data
+
+## In this case, the countries with minimal data (< 4 indicators ever calculated) are deleted.
+## These will be gap-filled later on if they are deleted now.
+d_gap_fill <- d_gap_fill %>%
+ filter(!(country %in% countries_no_data))This involves:
+d_calcs <- d_gap_fill %>%
+ group_by(country, year) %>%
+ summarize(score_wgi_scale = mean(value, na.rm=T),
+ NA_start = mean(NA_count_c_y), # initial mean number of NA across indicators, pre-gapfill
+ NA_post_gf_1 = mean(NA_count_post_gf1)) %>% # number of NA across indicators, post-gapfill across year gapfill within region/indicator
+ ungroup() scores_wgi_scale fall within the wgi range specified
+below:summary(d_calcs) # checking to make sure NA values make sense
+ ## v2022: NA_start max = 6, NA_post_gf_1 = 3
+
+## v2022: score_wgi_sclae min = -2.45, ma = 1.97... perfect
+wgi_range = c(-2.5, 2.5) # historic values have been between -2.5 and 2.5d_calcs <- d_calcs %>%
+ mutate(score = (score_wgi_scale - wgi_range[1]) / (wgi_range[2] - wgi_range[1])) %>%
+ ungroup(); head(d_calcs); summary(d_calcs)
+# converts the scores between 0 and 1, using min/max to scale the data
+
+# document gapfilling
+d_calcs <- d_calcs %>%
+ mutate(gap_fill = NA_start - NA_post_gf_1, # if there are values in NA_post_gf_1, it means these weren't gapfilled
+ gap_fill = ifelse(is.na(score), 0, gap_fill)) %>% # number of values that were gapfilled
+ select(-NA_start, -NA_post_gf_1)d_calcs[d_calcs$gap_fill>0, ] %>% View()
+d_calcs[d_calcs$country == "Niue", ] %>% View() # should have gap-fill values between 0-6
+d_calcs[d_calcs$country == "American Samoa",] %>% View() # should have gap-fill values between 0-6## save intermediate file of wgi scores pre-gapfilling (for OHI+ use)
+write.csv(d_calcs %>%
+ select(country, year, score_wgi_scale, score_ohi_scale = score),
+ here(paste0('globalprep/prs_res_wgi/v', scen_year, '/intermediate/wgi_combined_scores_by_country.csv')),
+ row.names = FALSE)## We report these regions at a greater spatial resolution:
+
+## Aruba is part of the Netherlands Antilles, but it is reported separately
+country_split_1 <- data.frame(country = "Netherlands Antilles", region = c('Bonaire', 'Curacao', 'Saba', 'Sint Maarten', 'Sint Eustatius'))
+country_split_2 <- data.frame(country = "Jersey, Channel Islands", region = c('Jersey', 'Guernsey'))
+country_split <- rbind(country_split_1, country_split_2)
+
+country_split_data <- country_split %>%
+ left_join(d_calcs) %>%
+ select(-country) %>%
+ rename(country = region)
+
+d_calcs <- d_calcs %>%
+ filter(!(country %in% c("Netherlands Antilles", "Jersey, Channel Islands"))) %>%
+ rbind(country_split_data) %>%
+ mutate(country = as.character(country))
+
+d_calcs$country[grep("Korea, Dem.", d_calcs$country)] <- "North Korea"
+
+d_calcs$country[grep("Turkiye", d_calcs$country)] <- "Turkey"
+# Maybe in future update package with country synonym list
+
+
+## Function to convert to OHI region ID
+d_calcs_rgn <- name_2_rgn(df_in = d_calcs,
+ fld_name='country',
+ flds_unique=c('year'))
+
+
+## Combine the duplicate regions (we report these at lower resolution)
+## In this case, we take the weighted average
+population_weights <- data.frame(country = c("Virgin Islands (U.S.)", "Puerto Rico",
+ "China", "Hong Kong SAR, China", "Macao SAR, China"),
+ population = c(106290, 3281538,
+ 1410929360, 7481800, 649342))
+# updated population values on 29 June 2022 (source: World Bank website, 2020 values)
+
+
+d_calcs_rgn <- d_calcs_rgn %>%
+ left_join(population_weights, by="country") %>%
+ mutate(population = ifelse(is.na(population), 1, population)) %>%
+ group_by(rgn_id, year) %>%
+ summarize(score = weighted.mean(score, population),
+ gapfill_within_rgn = weighted.mean(gap_fill, population)) %>%
+ ungroup() %>%
+ filter(rgn_id <= 250)
+
+summary(d_calcs_rgn)Assigning territorial region value to be the mean of parent country
+value and territorial regions with data (using same
+sov_id).
## data that describes territories of countries
+territory = rgn_master %>%
+ select(rgn_id = rgn_id_2013,
+ sov_id) %>%
+ group_by(rgn_id) %>% # remove duplicated countries from this rgn_id list
+ summarize(sov_id = mean(sov_id, na.rm=T)) %>% # duplicates always have the same sov_id (r2 value)
+ filter(rgn_id <= 250, rgn_id != 213)
+
+
+## expand to include all years of data
+territory <- data.frame(year=yr_start:yr_end) %>%
+ merge(territory, by=NULL)
+
+
+## assign territories the values of their sovereign country
+d_sovs = d_calcs_rgn %>%
+ full_join(territory, by = c('rgn_id', 'year')) %>%
+ group_by(sov_id, year) %>%
+ mutate(score_gf_territory = mean(score, na.rm=TRUE),
+ gapfill_within_rgn = mean(gapfill_within_rgn, na.rm=TRUE))%>%
+ filter(!is.na(gapfill_within_rgn)) %>%
+ ungroup()
+
+# filter(d_sovs, rgn_id %in% c(1,2,3,16) & year == 2018) %>% View() # check in console to make sure Australia and territories have the same score and sov_idDefine new data object from d_sovs which includes gapfill method and +gapfilled scores:
+d_gf2 <- d_sovs %>%
+ mutate(gapfill_territory = ifelse(is.na(score) & !is.na(score_gf_territory), "territory", "NA")) %>%
+ mutate(score = ifelse(is.na(score), score_gf_territory, score)) %>%
+ select(rgn_id, year, score, gapfill_within_rgn, gapfill_territory)
+
+summary(d_gf2)Add region names and clean the region data, and make sure we have all +the regions:
+# get region names
+regions <- rgn_master %>%
+ filter(rgn_typ == "eez") %>%
+ filter(rgn_id_2013 <= 250) %>% # > 250 are either FAO or a disputed region
+ filter(rgn_id_2013 != 213) %>% # 213 is Antarctica
+ select(rgn_id = rgn_id_2013, rgn_name=rgn_nam_2013) %>%
+ unique() %>%
+ arrange(rgn_id)
+
+d_gf2 <- regions %>%
+ left_join(d_gf2)head(d_sovs)
+summary(d_sovs)
+ ## v2021: score has 1029 NAs
+ ## v2022: 1127 NAs## check for NA values within "score" variable
+## if so, need to gapfill using UN geopolitical regions
+summary(d_gf2)
+d_gf2[is.na(d_gf2$score), ]
+ ## v2021: score has 0 NAs
+ ## v2022: 0 NAsComparing this year’s values against last year’s. These should be the +same unless there have been updates to WGI source data or a change to +methods. For this year (2019), there was a small change that affected a +few territorial regions. In the past, we used the sovereign country +value, but in this case, we averaged the sovereign country and the +available territorial values.
+Plot most recent shared year between last year and this year’s data, +and look for a relationship close to a 1:1 relationship. If data are +significantly off the line, look at the original (raw) data to +investigate.
+# Most recent shared year should be three years before prior to the assessment year (scen_year)
+recent_shared_year <- as.numeric(scen_year) - 3
+prev_assessment_year <- as.numeric(scen_year) - 1
+
+new_2019_scores <- d_gf2 %>% # rename variable with respective year for clarity
+ filter(year == recent_shared_year) %>%
+ select(rgn_id, score)
+
+old_2019_scores <- read.csv(here(paste0('globalprep/prs_res_wgi/v', prev_assessment_year, '/output/wgi_res.csv'))) %>% # rename variable
+ filter(year == recent_shared_year) %>%
+ select(rgn_id, old_score=resilience_score)
+
+score_compare_2019 <- old_2019_scores %>%
+ full_join(new_2019_scores)
+
+
+score_compare_plot <- ggplot(score_compare_2019, aes(x = old_score, y = score, text = rgn_id)) +
+ geom_point() +
+ geom_abline(slope=1, intercept = 0, col = "red")
+
+ggplotly(score_compare_plot) # looks good (no outliers), no real need for comparisonLook at top/bottom 10 regions to make sure these seem reasonable:
+## Top/Bottom 10 scorers:
+
+tmp <- d_gf2 %>%
+ filter(year==2019) %>%
+ arrange(score) %>%
+ select(rgn_id, score) %>%
+ left_join(regions)
+
+tmp[1:10, ] # makes sense.. mostly African countries or countries with authoritarian regimes
+tmp[211:220, ] # makes sense... mostly countries with democratic socialist policies or small islands. No surprise New Zealand is the best
+
+hist(tmp$score)Look at a summary to confirm scores are between 0 and 1, there are +220 regions, and there are no NAs (for this particular dataset):
+summary(d_gf2) # scores are between 0-1, gapfilled info present, no NAs
+length(unique(d_gf2$rgn_id)) # 220 regions
+c(min(d_gf2$score), max(d_gf2$score)) # checking for score min and max; making sure they're between 0-1Save gapfilling and data for this assessment year.
+tmp_data_res <- d_gf2 %>%
+ select(rgn_id, year, resilience_score = score)
+write.csv(tmp_data_res, here(paste0("globalprep/prs_res_wgi/v", scen_year, "/output/wgi_res.csv")), row.names=FALSE)
+
+tmp_data_prs <- d_gf2 %>%
+ mutate(score = 1 - score) %>%
+ select(rgn_id, year, pressure_score = score)
+write.csv(tmp_data_prs, here(paste0("globalprep/prs_res_wgi/v", scen_year, "/output/wgi_prs.csv")), row.names=FALSE)
+
+# gapfilled data
+tmp_gf <- d_gf2 %>%
+ select(rgn_id, year, gapfill_within_rgn, gapfill_territory) %>%
+ mutate(gapfill_within_rgn = ifelse(gapfill_within_rgn == 0, NA, gapfill_within_rgn)) %>%
+ mutate(gapfill_within_rgn = ifelse(!is.na(gapfill_within_rgn),
+ paste0("gapfill_within_rgn: ", gapfill_within_rgn), NA)) %>%
+ mutate(gapfill_territory = ifelse(gapfill_territory == "territory", "territory, mean of admin countries", NA)) %>%
+ mutate(method = paste(gapfill_within_rgn, gapfill_territory, sep="; ")) %>%
+ mutate(method = gsub("NA; ", "", method)) %>%
+ mutate(method = gsub("; NA", "", method)) %>%
+ mutate(method = gsub("NA", NA, method)) %>%
+ mutate(gapfilled = ifelse(is.na(method), 0, 1)) %>%
+ select(rgn_id, year, gapfilled, method)
+
+
+write.csv(tmp_gf, here(paste0("globalprep/prs_res_wgi/v", scen_year, "/output/wgi_gf.csv")), row.names=FALSE)Checking on outliers in 2022 assessment:
+## Do this after updating in global
+
+worldbank_wgi_from_wdi_api <- read_csv(paste0("globalprep/prs_res_wgi/v", scen_year, "/raw/worldbank_wgi_from_wdi_api_2022-10-03.csv"))
+
+worldbank_2018 <- worldbank_wgi_from_wdi_api %>%
+ filter(year == 2018) %>%
+ filter(country == "New Zealand")
+
+# Compare to data from 2021:
+worldbank_wgi_from_wdi_api_2021 <- read_csv(paste0("globalprep/prs_res_wgi/v", prev_assessment_year, "/raw/worldbank_wgi_from_wdi_api_2021-04-29.csv"))
+
+worldbank_2018_old <- worldbank_wgi_from_wdi_api_2021 %>%
+ filter(year == 2018) %>%
+ filter(country == "New Zealand")