Replication script for a chapter on the economic impacts of Malaria.

malara_elimination.R

# Country map ---

library("sf")
library("dplyr")
library("tmap")

data("World")
World_o <- World
# Merge Somalia and Somaliland
World <- World |>
  mutate(iso_a3 = as.character(iso_a3), name = as.character(name)) |>
  mutate(iso_a3 = ifelse(iso_a3 == "SOL", "SOM", iso_a3)) |>
  mutate(name = ifelse(name == "Somaliland", "Somalia", name)) |>
  group_by(iso_a3) |> summarise() |>
  left_join(World |> st_drop_geometry(), by = "iso_a3")

d_eli <- readr::read_csv("data/elim_countries.csv")
d_inc <- readr::read_csv("data/incidence_countries.csv", skip = 3) |>
  select(iso_a3 = "Country Code", as.character(2000:2018))
d_inc <- d_inc[!apply(d_inc[, -1], 1, \(x) all(is.na(x))), ]
d_inc$incidence <- cut(d_inc[["2018"]], c(0, 1, 10, 100, 250, Inf))

data <- World |>
  left_join(d_inc, by = "iso_a3") |>
  left_join(d_eli, by = "name") |>
  mutate(value = ifelse(is.na(status), as.character(incidence), status)) |>
  mutate(free = ifelse(status == "greater_one", "under 3 years", "3 to 20 years"))
# data$free <- addNA(data$free)
# levels(data$free)[3] <- "more than 20 years" # very dirty

map <- data |>
  rename(`Malaria free` = free) |>
  st_transform(crs = "+proj=moll +x_0=0 +y_0=0 +lat_0=0 +lon_0=0")

t <- tm_shape(map, bbox = c(-12000000, -5500000, 16000000, 7500000)) +
  tm_fill("incidence", palette =
    rev(c("#3B3B3B", "#575757", "#ABABAB", "#DADADA", "#EBEBEB")),
    # rev(c("#7D0025", "#b30000", "#fc8d59", "#fdd49e", "#fee8c8")), # Coloured
    labels = c("0 to 0.1", "0.1 to 1", "1 to 10", "10 to 25", "25 to 100"),
    title = "Incidence (%)", showNA = FALSE, colorNA = NULL) +
    tm_borders(col = "#bbbbbb") +
  tm_shape(map) +
  # tm_graticules(labels.size = 1.2, lwd = 0.25, n.y = 1, n.x = 2) +
  tm_symbols(shape = "Malaria free", shapes = c(15, 17), col = c("#333333"),
    shapeNA = NA, shape.showNA = FALSE,
    # palette = c("#9ecae1", "#deebf7", "#f4f4f4"),
    # shapes.legend = "Malaria free",
    labels = c("<3 years", "3 to 20 years")) +
    tm_layout(frame = TRUE, fontfamily = "Helvetica",
    title = "Malaria status",
    title.position = c("center", "top"), title.size = 2.4,
    title.fontface = "bold",
    outer.bg.color = "transparent", bg.color = "transparent",
    inner.margins = c(0.01, 0.01, 0.01, 0.01),
    legend.position = c("left", "bottom"), legend.frame = FALSE,
    legend.text.size = 1.2, legend.title.size = 1.8)
t

tmap_save(t, "incidence_map.pdf", dpi = 140,
  height = 950, width = 2000, bg = "transparent")
# tmap_save(t, "incidence_map.png", dpi = 140,
#   height = 950, width = 2000, bg = "transparent")


# Subnational map ---

# Add Venezuela's poverty data
ven <- read_sf("data/gadm36_gpkg/gadm36.gpkg")
ven <- ven |> dplyr::filter(GID_0 == "VEN") |>
  dplyr::group_by(NAME_1) |> dplyr::summarise()
ven <- ven |>
  dplyr::left_join(readr::read_csv("data/poverty_ven.csv"),
    by = c("NAME_1" = "state"))
ven <- st_transform(ven, crs = "+proj=moll +x_0=0 +y_0=0 +lat_0=0 +lon_0=0")
saveRDS(ven, "data/poverty_ven.rds")

library("stars")
library("sf")
library("tmap")

data("World")
World <- st_transform(World, crs = "+proj=moll +x_0=0 +y_0=0 +lat_0=0 +lon_0=0")

r1 <- stars::read_stars("data/Pf_Incidence/Raster Data/Pf_incidence_rate_rmean/2020_GBD2019_Global_Pf_Incidence_Rate_2019.tif")
r2 <- stars::read_stars("data/Pv_Incidence/Raster Data/Pv_incidence_rate_rmean/2020_GBD2019_Global_Pv_Incidence_2019.tif")
s <- sf::st_read("data/GSAP2/GSAP2.shp")
ven <- readRDS("data/poverty_ven.rds")

s <- st_transform(s, crs = "+proj=moll +x_0=0 +y_0=0 +lat_0=0 +lon_0=0")
r1 <- st_transform(r1, st_crs(s))
r2 <- st_transform(r2, st_crs(s))

#t <- s |> dplyr::mutate(value = as.numeric(GSAP2_poor)) |>
#  tm_shape(bbox = c(-12000000, -5500000, 16000000, 7500000)) +
#  tm_fill("value", palette = "BuPu", style = "fixed",
#          title = "Rate (%)",
#          breaks = c(0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 100)) +
#  tm_shape(ven |> dplyr::mutate(value = poor_extr)) +
#  tm_fill("value", palette = "BuPu", style = "fixed",
#          legend.show = FALSE,
#          breaks = c(0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 100)) +
#  tm_shape(s) +
#  tm_borders(col = "#dddddd", alpha = 0.5) +
#  tm_shape(World) + tm_borders(col = "#bbbbbb") +
#  tm_layout(frame = TRUE, fontfamily = "Helvetica",
#            title = "Extreme poverty",
#            title.position = c("center", "top"), title.size = 2.4,
#            title.fontface = "bold",
#            outer.bg.color = "transparent", bg.color = "transparent",
#            inner.margins = c(0.01, 0.01, 0.01, 0.01),
#            legend.position = c("left", "bottom"), legend.frame = FALSE,
#            legend.text.size = 1.2, legend.title.size = 1.8)

#tmap_save(t, "poverty.pdf", dpi = 140,
#          height = 950, width = 2000, bg = "transparent")
#tmap_save(t, "poverty.png", dpi = 140,
#          height = 950, width = 2000, bg = "transparent")


t <- tm_shape(r1 * 100, raster.downsample = FALSE,
              bbox = c(-12000000, -5500000, 16000000, 7500000)) +
  # tm_graticules(labels.size = 1.2, lwd = 0.25, n.y = 1, n.x = 2) +
  tm_raster(palette = "OrRd", style = "fixed",
            title = "Incidence (%)",
            breaks = c(0, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 1) * 100) +
  tm_shape(World) + tm_borders(col = "#bbbbbb") +
  tm_layout(frame = TRUE, fontfamily = "Helvetica",
            title = "P. falciparum",
            title.position = c("center", "top"), title.size = 2.4,
            title.fontface = "bold",
            outer.bg.color = "transparent", bg.color = "transparent",
            inner.margins = c(0.01, 0.01, 0.01, 0.01),
            legend.position = c("left", "bottom"), legend.frame = FALSE,
            legend.text.size = 1.2, legend.title.size = 1.8)

#tmap_save(t, "incidence_falc.pdf", dpi = 140,
#          height = 950, width = 2000, bg = "transparent")
tmap_save(t, "incidence_falc.png", dpi = 140,
          height = 950, width = 2000, bg = "transparent")


t <- tm_shape(r2 * 100, raster.downsample = FALSE,
              bbox = c(-12000000, -5500000, 16000000, 7500000)) +
  # tm_graticules(labels.size = 1.2, lwd = 0.25, n.y = 1, n.x = 2) +
  tm_raster(palette = "OrRd", style = "fixed",
            title = "Incidence (%)",
            breaks = c(0, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 1) * 100) +
  tm_shape(World) + tm_borders(col = "#bbbbbb") +
  tm_layout(frame = TRUE, fontfamily = "Helvetica",
            title = "P. vivax",
            title.position = c("center", "top"), title.size = 2.4,
            title.fontface = "bold",
            outer.bg.color = "transparent", bg.color = "transparent",
            inner.margins = c(0.01, 0.01, 0.01, 0.01),
            legend.position = c("left", "bottom"), legend.frame = FALSE,
            legend.text.size = 1.2, legend.title.size = 1.8)

#tmap_save(t, "incidence_vivax", dpi = 140,
#          height = 950, width = 2000, bg = "transparent")
tmap_save(t, "incidence_vivax.png", dpi = 140,
          height = 950, width = 2000, bg = "transparent")


# Country map ---

library("sf")
library("dplyr")
library("tmap")

data("World")
World_o <- World
# Merge Somalia and Somaliland
World <- World %>%
  mutate(iso_a3 = as.character(iso_a3), name = as.character(name)) %>%
  mutate(iso_a3 = ifelse(iso_a3 == "SOL", "SOM", iso_a3)) %>%
  mutate(name = ifelse(name == "Somaliland", "Somalia", name)) %>%
  group_by(iso_a3) %>% summarise() %>%
  left_join(World %>% st_drop_geometry(), by = "iso_a3")

d_code <- readr::read_csv("data/code_spending.csv")
d_spending <- readr::read_csv("data/data_spending.csv")

d_spending <- d_spending %>%
  filter(level == "Country", year == 2016) %>%
  select(iso3_rc, location_name, year, the_total_mean, the_per_case_mean) %>%
  rename(iso3c = iso3_rc, country = location_name, year = year,
         tot_mean = the_total_mean, pc_mean = the_per_case_mean) %>%
  left_join(d_code %>% select(iso3c, status), by = "iso3c") %>%
  mutate(pc_mean_e = ifelse(pc_mean * (status == "e") == 0, NA, pc_mean * (status == "e")),
         pc_mean_c = ifelse(pc_mean * (status == "c") == 0, NA, pc_mean * (status == "c")),
         pc_mean_f = ifelse(pc_mean * (status == "f") == 0, NA, pc_mean * (status == "f")),
         pc_mean_e_bins = case_when(
           pc_mean_e <= 10 ~ "0 to 10",
           pc_mean_e > 10  & pc_mean_e <= 50 ~ "10 to 50",
           pc_mean_e > 50  & pc_mean_e <= 250 ~ "50 to 250",
           pc_mean_e > 250  & pc_mean_e <= 2500 ~ "250 to 2,500",
           pc_mean_e > 2500 ~ "2,500+"),
         pc_mean_c_bins = case_when(
           pc_mean_c <= 10 ~ "0 to 10",
           pc_mean_c > 10  & pc_mean_c <= 50 ~ "10 to 50",
           pc_mean_c > 50  & pc_mean_c <= 250 ~ "50 to 250",
           pc_mean_c > 250  & pc_mean_c <= 2500 ~ "250 to 2,500",
           pc_mean_c > 2500 ~ "2,500+"))

d_spending$pc_mean_e_bins <- factor(d_spending$pc_mean_e_bins,
                                    levels = c("0 to 10", "10 to 50",
                                               "50 to 250", "250 to 2,500",
                                               "2,500+"))

d_spending$pc_mean_c_bins <- factor(d_spending$pc_mean_c_bins,
                                    levels = c("0 to 10", "10 to 50",
                                               "50 to 250", "250 to 2,500",
                                               "2,500+"))
data <- World %>%
  left_join(d_spending, by = c("iso_a3" = "iso3c"))

map <- data %>%
  mutate(Status = ifelse(!is.na(pc_mean_c_bins), "Controlling",
    ifelse(!is.na(pc_mean_e_bins), "Eliminating", NA))) %>%
  st_transform(crs = "+proj=moll +x_0=0 +y_0=0 +lat_0=0 +lon_0=0")

t <- tm_shape(map, bbox = c(-12000000, -5500000, 16000000, 7500000)) +
  tm_fill("pc_mean_c_bins", style = "fixed",
          palette = "Greys",
          labels = c("0 to 10", "10 to 50", "50 to 250", "250 to 2,500", "2,500+"),
          legend.show = FALSE, title = "Controlling",
          showNA = FALSE, colorNA = NULL, shapes = 4) +
  tm_shape(map) +
  tm_fill("pc_mean_e_bins", style = "fixed",
          palette = "Greys",
          labels = c("0 to 10", "10 to 50", "50 to 250", "250 to 2,500", "2,500+"),
          title = "Spending",
          showNA = FALSE, colorNA = NULL) +
  tm_borders(col = "#bbbbbb") +
  tm_symbols(shape = "Status", shapes = c(21, 24), size = .5,
    col = c("#333333"), border.col = "#a4a4a4",
    shapeNA = NA, shape.showNA = FALSE) +
  tm_layout(frame = TRUE, fontfamily = "Helvetica",
            title = "Spending per malaria case",
            title.position = c("center", "top"), title.size = 2.4,
            title.fontface = "bold",
            outer.bg.color = "transparent", bg.color = "transparent",
            inner.margins = c(0.01, 0.01, 0.01, 0.01),
            legend.position = c("left", "bottom"), legend.frame = FALSE,
            legend.text.size = 1.2, legend.title.size = 1.8)

tmap_save(t, "spending_combined.pdf", dpi = 140,
  height = 950, width = 2000, bg = "transparent")