library(tidyverse)
library(scales)
library(sf)
library(patchwork)
library(countrycode)
library(readxl)
library(validate)
library(pander)
library(here)

source(here("lib", "graphic-functions.R"))
# qwraps2::lazyload_cache_dir("figures_cache/html")

options(dplyr.summarise.inform = FALSE)

CIVICUS map

# Load CIVICUS data
# Downloaded from HTML page at https://monitor.civicus.org/
civicus <- read_csv(here("data", "civicus_monitor_2017.csv"), na = "Null") %>%
  mutate(Rating = factor(Rating, levels = c("Open", "Narrowed", "Obstructed", 
                                            "Repressed", "Closed"), 
                         ordered = TRUE),
         iso3 = countrycode(Country, "country.name", "iso3c"))

# Load global shapefiles
# Get Admin 0 file from http://www.naturalearthdata.com/downloads/110m-cultural-vectors/
world_shapes <- st_read(here("data", "ne_110m_admin_0_countries",
                             "ne_110m_admin_0_countries.shp"),
                        quiet = TRUE) %>% 
  # Ignore Antarctica
  filter(ISO_A3 != "ATA")

# Use the Robinson map projection
projection = "ESRI:54030"

# Join CIVICUS data with map data
map_with_civicus <- world_shapes %>% 
  # Fix some Natural Earth ISO weirdness
  mutate(ISO_A3 = ifelse(ISO_A3 == "-99", as.character(ISO_A3_EH), as.character(ISO_A3))) %>% 
  mutate(ISO_A3 = case_when(
    .$ISO_A3 == "GRL" ~ "DNK",
    .$NAME == "Norway" ~ "NOR",
    TRUE ~ ISO_A3
  )) %>% 
  left_join(civicus, by = c("ISO_A3" = "iso3"))

plot_civicus_map <- ggplot() +
  geom_sf(data = map_with_civicus, aes(fill = Rating), size = 0.15, color = "black") +
  coord_sf(crs = st_crs(projection), datum = NA) +
  scale_fill_manual(values = c("grey90", "grey70", "grey45",
                               "grey20", "black"),
                    na.translate = FALSE, name = "Civic space") +
  theme_ngo_map(9)
plot_civicus_map

# Save the map
ggsave(plot_civicus_map, filename = here("output", "civicus_map.pdf"), 
       width = 5.5, height = 3, units = "in", device = cairo_pdf)

ggsave(plot_civicus_map, filename = here("output", "civicus_map.png"), 
       width = 5.5, height = 3, units = "in", dpi = 300, type = "cairo")

# ggsave(plot_civicus_map, filename = here("output", "civicus_map.tiff"), 
#        width = 5.5, height = 3, units = "in", type = "cairo", dpi = 600)

Donors and restrictions

predicted_oda <- read_csv(here("data", "predicted-donors",
                               "predicted_h1_barriers_individual.csv"))

predicted_contention <- read_csv(here("data", "predicted-donors",
                                      "predicted_h2_barriers_individual.csv"))
df_oda_predict <- predicted_oda %>% 
  gather(barrier, barrier_value, advocacy_within) %>%
  group_by(cowcode, barrier_value, barrier) %>%
  summarise_at(vars(predicted), list(~mean(.))) %>% 
  mutate(fake_facet_title = "Barriers to advocacy and aid")

df_oda_predict_mean <- df_oda_predict %>%
  group_by(barrier, barrier_value) %>%
  summarise(predicted = mean(predicted)) %>%
  mutate(cowcode = 1)  # Fake country so the line plots

plot_oda_predict <- ggplot(df_oda_predict, 
       aes(x = barrier_value, y = expm1(predicted), group = cowcode)) + 
  geom_vline(xintercept = 0, colour = "black", size = 0.5) + 
  geom_smooth(method = "lm", size = 0.1, alpha = 0.1, 
              colour = "grey30") +
  geom_smooth(data = df_oda_predict_mean, size = 1.5, 
              method = "lm", colour = "black") +
  labs(x = "Change from average number of anti-advocacy barriers",
       y = "Predicted ODA in the following year") +
  scale_y_continuous(labels = dollar) +
  theme_ngo() +
  facet_wrap(~ fake_facet_title)

# Inverse logit, with the ability to account for adjustments
# via http://stackoverflow.com/a/23845527/120898
inv_logit <- function(f, a) {
  a <- (1 - 2 * a)
  (a * (1 + exp(f)) + (exp(f) - 1)) / (2 * a * (1 + exp(f)))
}

df_contention_predict <- predicted_contention %>% 
  gather(barrier, barrier_value, advocacy_within) %>%
  group_by(cowcode, barrier_value, barrier) %>%
  summarise_at(vars(predicted), list(~mean(.))) %>% 
  mutate(predicted_fixed = inv_logit(predicted, a = 0.001)) %>% 
  mutate(fake_facet_title = "Barriers to advocacy\nand contentiousness of aid")

df_contention_predict_mean <- df_contention_predict %>%
  group_by(barrier, barrier_value) %>%
  summarise(predicted = mean(predicted)) %>%
  mutate(cowcode = 1,  # Fake country so the line plots
         predicted_fixed = inv_logit(predicted, a = 0.001))

plot_contention_predict <- ggplot(df_contention_predict, 
       aes(x = barrier_value, y = predicted_fixed, group = cowcode)) + 
  geom_vline(xintercept = 0, colour = "black", size = 0.5) + 
  geom_smooth(method = "lm", size = 0.1, alpha = 0.1, 
              colour = "grey30") +
  geom_smooth(data = df_contention_predict_mean, size = 1.5, 
              method = "lm", colour = "black") +
  labs(x = "Change from average number of anti-advocacy barriers",
       y = "Predicted proportion of\ncontentious aid in the following year") +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  theme_ngo() +
  facet_wrap(~ fake_facet_title)

plot_predict_both <- 
  (plot_oda_predict + plot_contention_predict + plot_layout(ncol = 2)) &
  theme(axis.title.x = element_text(size = rel(0.8)),
        axis.title.y = element_text(margin = margin(r = 3), size = rel(0.8)))

plot_predict_both

ggsave(plot_predict_both, filename = here("output", "fig-predict-both.pdf"), 
       width = 5.5, height = 2.5, units = "in", device = cairo_pdf)

ggsave(plot_predict_both, filename = here("output", "fig-predict-both.png"), 
       width = 5.5, height = 2.5, units = "in", dpi = 300, type = "cairo")

# ggsave(plot_predict_both, filename = here("output", "fig-predict-both.tiff"), 
#        width = 5.5, height = 2.5, units = "in", type = "cairo", dpi = 600)