4_SchoolReforms_Add_EMDAT_disaster_data.R

###############################################################
###############################################################
#######                                                 ####### 
#######                   PROJECT:                      ####### 
#######               School reforms                    ####### 
#######                                                 ####### 
#######           CODE: Add climate data                ####### 
###############################################################
###############################################################
rm(list =ls())

library(sf)
library(raster)
library(viridis)
library(rgdal)
library(rworldmap)
library(tidyverse)
library(lubridate)
library(dplyr)
library(ggplot2)
#install.packages("here")
library(here)


here::here()    #shows where the project root is.
#here::dr_here() #tells you why 

#####################################################################
### 1. Prepare the DHS data
#####################################################################

data <- read.csv("D:/Anna/Dropbox/Projects/2024_San Diego/School reforms/Data/data_for_analysis.csv")
data <- data %>% dplyr::select(-c(X))
cntr <- unique(data$CountryName)


# list of countries included in the analysis
cntr_control = c("Burkina Faso", "Congo Democratic Republic", "Lesotho", "Liberia", "Zimbabwe")
cntr_exposed = c("Ethiopia", "Kenya", "Mali", "Senegal", "Tanzania", "Uganda")
cntr <- c(cntr_control, cntr_exposed)

data <- data %>% 
  #convert Ethiopian calendar months to numeric
  mutate(intMo = ifelse(intMo=="magabit", 7, intMo)) %>% 
  mutate(intMo = ifelse(intMo=="miyazya", 8, intMo)) %>% 
  mutate(intMo = ifelse(intMo=="ginbot", 9, intMo)) %>% 
  mutate(intMo = ifelse(intMo=="sene", 10, intMo)) %>% 
  #convert other months to numeric
  mutate(intMo = ifelse(intMo=="june", 6, intMo)) %>% 
  mutate(intMo = ifelse(intMo=="july", 7, intMo)) %>% 
  mutate(intMo = ifelse(intMo=="august", 8, intMo)) %>% 
  mutate(intMo = ifelse(intMo=="september", 9, intMo)) %>% 
  mutate(intMo = ifelse(intMo=="october", 10, intMo)) %>% 
  mutate(intMo = ifelse(intMo=="november", 11, intMo)) %>% 
  mutate(intMo = as.numeric(intMo)) %>% 
  mutate(intYr = as.numeric(intYr)) %>% 
  # convert dates from Ethiopian to Gregorian calendar
  mutate(intCMC   = ifelse(CountryName == "Ethiopia", intCMC+92, intCMC)) %>% 
  mutate(birthCMC = ifelse(CountryName == "Ethiopia", birthCMC+92, birthCMC)) %>% 
  # using cmc to determine year and month in DHS: https://dhsprogram.com/data/Guide-to-DHS-Statistics/index.cfm
  mutate(intYr_new   = as.integer((intCMC - 1)/12)+1900) %>% 
  mutate(intMo_new   = intCMC - ((intYr_new - 1900)*12)) %>% 
  mutate(birthYr_new = as.integer((birthCMC - 1)/12)+1900) %>% 
  mutate(birthMo_new = birthCMC - ((birthYr_new - 1900)*12)) %>% 
  mutate(intYr       = ifelse(CountryName == "Ethiopia", intYr_new, intYr)) %>% 
  mutate(intMo       = ifelse(CountryName == "Ethiopia", intMo_new, intMo)) %>% 
  mutate(birthYr     = ifelse(CountryName == "Ethiopia", birthYr_new, birthYr)) %>% 
  mutate(birthMo     = ifelse(CountryName == "Ethiopia", birthMo_new, birthMo)) %>%     
  dplyr::select(-c(intYr_new, intMo_new, birthYr_new, birthMo_new))


# Interview year missing for Congo DRC - correct
data <- data  %>% 
  mutate(intYr_new   = as.integer((intCMC - 1)/12)+1900) %>% 
  mutate(intYr       = ifelse(CountryName == "Congo Democratic Republic", intYr_new, intYr)) %>%     
  dplyr::select(-c(intYr_new)) %>% 
  filter(CountryName %in% cntr)


psu <- data %>% 
  dplyr::select(CountryName, SurveyId, psu, LATNUM, LONGNUM) %>% 
  unique()

unique(psu$CountryName)

#####################################################################
### 2. Prepare the EMDAT data
#####################################################################

load("D:/Anna/Dropbox/Data/pend-gdis-1960-2018-disasterlocations-rdata/pend-gdis-1960-2018-disasterlocations.Rdata")
str(GDIS_disasterlocations)
class(GDIS_disasterlocations)

sort(unique(GDIS_disasterlocations$country))

GDIS <- GDIS_disasterlocations %>% 
  mutate(country = ifelse(country=="Democratic Republic Of The Congo", "Congo Democratic Republic", country)) %>% 
  filter(country %in% cntr)
sort(unique(GDIS$country))

emdat <- read.csv("D:/Anna/Dropbox/Data/EM-DAT/public_emdat_custom_request_2025-10-20_568efae0-5382-417b-b95a-5ae410862b83.csv")
unique(emdat$Country)
emdat <- emdat %>% 
  mutate(Country = ifelse(Country=="Democratic Republic of the Congo", "Congo Democratic Republic", Country)) %>% 
  mutate(Country = ifelse(Country=="United Republic of Tanzania", "Tanzania", Country)) %>% 
  filter(Country %in% cntr)
sort(unique(emdat$Country))
names(emdat)[1] <- "Dis.No"
unique(emdat$Disaster.Subgroup)

emdat <- emdat %>% 
  #filter(Disaster.Group=="Natural") %>% 
  filter(Disaster.Subgroup == "Hydrological" | Disaster.Subgroup == "Meteorological" | Disaster.Subgroup == "Climatological") %>% 
  dplyr::select(c(Dis.No, Country, Subregion, Region, Location, Disaster.Group, Disaster.Subgroup, Disaster.Type, Disaster.Subtype, Event.Name, Location, Start.Year, Start.Month, Start.Day, End.Year, End.Month, End.Day, Total.Affected)) %>% 
  mutate(Dis.No = str_sub(Dis.No, end = -5)) 

emdat <- emdat %>% 
  filter(Start.Year >= 1993 & Start.Year <= 2018) %>% 
  filter(Dis.No %in% GDIS$disasterno) 

GDIS <- GDIS %>% 
  filter(disasterno %in% emdat$Dis.No)

unique(GDIS$disastertype)


#####################################################################
### 3. Map the EMDT data at sub-national level
#####################################################################
library(geodata)
library(sf)
library(dplyr)
library(purrr)

# define countries
cntr_codes <- c("BFA", "COD", "LSO", "LBR", "ZWE",
               "ETH", "KEN", "MLI", "SEN", "TZA", "UGA")

# create directory if missing
if(!dir.exists("data")) dir.create("data")

# safe wrapper for gadm() that downloads if missing
get_admin2 <- function(iso) {
  message("Retrieving admin-2 boundaries for: ", iso)
  tryCatch({
    # explicitly specify version = "4.1" (latest) and download = TRUE
    g <- geodata::gadm(country = iso, level = 2, path = "data", version = "4.1")
    g_sf <- st_as_sf(g)
    g_sf$ISO3 <- iso
    g_sf
  }, error = function(e) {
    message("❌ Failed to retrieve admin-2 for ", iso, ": ", e$message)
    NULL
  })
}




# loop over all and bind
adm2 <- purrr::map_dfr(cntr_codes, get_admin2)

# transform CRS
adm2 <- st_transform(adm2, 4326)

unique(adm2$COUNTRY)

# Admin level 1 for Lesotho
adm1 <- geodata::gadm(country = "LSO", level = 1, path = "data/") |>
  st_as_sf() |>
  st_transform(4326)


plot(adm2)

adm2 <- bind_rows(adm2, adm1)
adm2 <- adm2 %>% 
  dplyr::select(GID_0, COUNTRY, GID_1, NAME_1, NAME_2, geometry) 

library(sf)
library(dplyr)

# Make sure CRS matches
adm2_sf <- st_make_valid(adm2)
GDIS_sf <- st_make_valid(GDIS)

adm2_sf <- st_transform(adm2_sf, 4326)      # WGS84
GDIS_sf <- st_transform(GDIS_sf, 4326)

adm2_sf <- adm2_sf %>% 
  # change to NAME_1 for Lesotho (no admin 2 available)
  mutate(NAME_2 = ifelse(is.na(NAME_2), NAME_1, NAME_2)) %>% 
  mutate(COUNTRY = ifelse(COUNTRY=="Democratic Republic of the Congo", "Congo Democratic Republic", COUNTRY))
  
# Overlay disasters onto admin-2
# Spatial join: for each admin-2 polygon, find disasters that intersect it
adm2_GDIS_hits <- st_join(
  adm2_sf,
  GDIS_sf %>% dplyr::select(disasterno, country, disastertype, geometry),
  join = st_intersects,
  left = TRUE
)


# Count unique disasters per admin-2
disaster_counts_admin2 <- adm2_GDIS_hits %>%
  st_drop_geometry() %>%                      # we just want attributes for counting
  group_by(NAME_2, country) %>%
  summarise(
    n_disasters = n_distinct(disasterno)      # number of unique disasters affecting this admin-2
  ) %>%
  ungroup()


# Join the counts back to the admin-3 polygons for mapping
adm2_with_counts <- adm2_sf %>%
  left_join(disaster_counts_admin2,
            by = c("COUNTRY"="country", "NAME_2"="NAME_2")) %>%
  mutate(n_disasters = ifelse(is.na(n_disasters), 0, n_disasters))

# Plot the map
# install.packages("rnaturalearth")
# install.packages("rnaturalearthdata")
library(rnaturalearth)
library(sf)
library(ggplot2)
library(viridis)


# Get background Africa map
africa_countries <- ne_countries(continent = "Africa", returnclass = "sf")

africa_countries_sub <- africa_countries %>% 
  filter(sovereignt %in% c(cntr, "Democratic Republic of the Congo", "United Republic of Tanzania"))

# Plot with Africa background first, then your admin2 polygons
map_1 <- ggplot() +
  # background layer in light gray
  geom_sf(data = africa_countries, fill = "gray70", color = NA, size = 0.01) + #color = "gray95"
  # admin2 data with disaster counts
  geom_sf(data = adm2_with_counts, aes(fill = n_disasters), color = NA, size = NA) +
  # color scale
  #scale_fill_viridis_c(
  #  option = "plasma",
  #  direction = 1,
  #  name = "Disasters") +
  scale_fill_distiller(
    palette = "Reds",
    direction = 1,
    name = "Disasters"
  )+ 
  #geom_sf(data = africa_countries_sub, fill = NA, color = "darkgray", size = 0.2) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    legend.position = "right") +
  labs(
    title = "Total Recorded Climate Disasters by Admin-2",
    subtitle = "1993–2018",
    x = NULL, y = NULL)


# Save as PNG
ggsave("disaster_map.png", plot = map_1,
       width = 10, height = 8, dpi = 300)

# Save as SVG
ggsave("disaster_map.svg", plot = map_1,
       width = 10, height = 8, device = "svg")



### Plot DHS location:
#srvy_list <- data %>% 
#  filter(intYr>=1992 & intYr<=2018) %>% 
#  group_by(CountryName, SurveyId) %>% summarise(count=n())
psu <- data %>% 
  filter(intYr>=1993 & intYr<=2018) %>% 
  dplyr::select(CountryName, SurveyId, psu, LATNUM, LONGNUM) %>% 
  unique()

pts_sf <- st_as_sf(psu, coords = c("LONGNUM","LATNUM"), crs = 4326, remove = FALSE)


# Plot with Africa background first, then the PSUs
map_2 <- ggplot() +
  # background layer in light gray
  geom_sf(data = africa_countries, fill = "gray70", color = "gray95", size = 0.01) + #color = "gray95"
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    legend.position = "right") +
  geom_point(
    data = pts_sf,
    aes(x = LONGNUM, y = LATNUM),
    size = 0.5, alpha = 0.9, na.rm = TRUE) +
  labs(x = NULL, y = NULL)

map_2

# Save as PNG
ggsave("PSUs_map.png", plot = map_2,
       width = 10, height = 8, dpi = 300)



### Plor disasters by type
unique(adm2_GDIS_hits$disastertype)

# Count unique disasters per admin-2 by type of disaster
disaster_counts_by_type_admin2 <- adm2_GDIS_hits %>%
  st_drop_geometry() %>%                      # we just want attributes for counting
  group_by(NAME_2, country, disastertype) %>%
  summarise(
    n_disasters = n_distinct(disasterno)      # number of unique disasters affecting this admin-2
  ) %>%
  ungroup()


# Join the counts back to the admin-3 polygons for mapping
adm2_with_counts_by_type <- adm2_sf %>%
  left_join(disaster_counts_by_type_admin2,
            by = c("COUNTRY"="country", "NAME_2"="NAME_2")) %>%
  mutate(n_disasters = ifelse(is.na(n_disasters), 0, n_disasters)) %>% 
  filter(!is.na(disastertype))

# Plot with Africa background first, then your admin2 polygons
map_3 <- ggplot() +
  # background layer in light gray
  geom_sf(data = africa_countries,
          fill = "gray70",
          color = "gray95",
          size = 0.01) +
  # admin-2 polygons with disaster counts
  geom_sf(data = adm2_with_counts_by_type,
          aes(fill = n_disasters),
          color = NA,
          size = 0.05) +
  # facet by disaster type — each facet gets its own title
  facet_wrap(~ str_to_title(disastertype),
             ncol = 2,
             labeller = label_value) +
  # color scale
  scale_fill_distiller(
    palette = "Reds",
    direction = 1,
    name = "Disasters"
  ) + 
  # overall layout and style
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    legend.position = "right",
    strip.background = element_rect(fill =NA, color = NA),
    strip.text = element_text(face = "bold", size = 12),
    plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
    plot.subtitle = element_text(size = 11, hjust = 0.5)
  ) +
  
  labs(
    #title = "Recorded Climate Disasters by Type and Admin-2 Region",
    #subtitle = "Across African countries, 1992–2018",
    x = NULL,
    y = NULL
  )

map_3

# Save as PNG
ggsave("disaster_by_type_map.png", plot = map_3,
       width = 10, height = 10, dpi = 300)


#####################################################################
### 4. Intersect the PSUs with the EMDAT data and generate exposures
#####################################################################

# Make points an sf object (WGS84). If you have a data.frame with lon/lat:
pts_sf <- st_as_sf(psu, coords = c("LONGNUM","LATNUM"), crs = 4326, remove = FALSE)
# Ensure polygons are sf and fix old-style CRS warnings
polys_sf <- GDIS  # your polygons
if (is.na(st_crs(polys_sf))) {
  # If GDIS came with an old proj4string but is actually WGS84, set it explicitly:
  st_crs(polys_sf) <- 4326
}
# Put both layers in the same CRS (use a projected CRS if you'll do distances)
polys_sf <- st_transform(polys_sf, st_crs(pts_sf))
# Fix invalid geometries to avoid join/intersection errors
polys_sf <- st_make_valid(polys_sf)


library(sf)
library(dplyr)
library(units)
# Spatial join: add polygon attributes to each point
# Use st_within for strict containment; st_intersects is more permissive.
#pts_in_polys <- st_join(pts_sf, polys_sf, join = st_within)
# distance-based predicate: within 10 km of polygon boundary (or inside)
within_10km <- function(x, y) st_is_within_distance(x, y, dist = set_units(10, km))
# Join polygon attributes to points if within 10 km (or inside)
pts_in_polys <- st_join(pts_sf, polys_sf, join = within_10km, left = TRUE)


# Now `pts_in_polys` is an sf POINT layer with polygon fields (e.g., country/adm1/.)
pts_in_polys <- pts_in_polys %>% st_drop_geometry()
pts_in_polys <- pts_in_polys %>% dplyr::select(-centroid)
pts_in_polys <-  pts_in_polys %>% 
  left_join(emdat, by=c("country"="Country", "disasterno" = "Dis.No")) %>% 
  dplyr::select(CountryName, SurveyId, psu, disastertype, disasterno, 
                Start.Year, Start.Month, End.Year, End.Month, Total.Affected)

#If end date is missing, use start date
#If start date is missing, drop
pts_in_polys <- pts_in_polys %>% 
  mutate(End.Month = ifelse(is.na(End.Month), Start.Month, End.Month))
  
pts_in_polys <- pts_in_polys %>% 
  filter(!is.na(Start.Year) & !is.na(Start.Month))

save(pts_in_polys, file = "Data/Emdat_data_extract.RData")

rm(list = setdiff(ls(), c("pts_in_polys", "data")))

#load("Data/Emdat_data_extract.RData")

### Create a data frame with the full time period and the disaster information
disas <- pts_in_polys %>% 
  dplyr::select(CountryName, SurveyId, psu, Start.Year, Start.Month, End.Year, End.Month)

# --- config: analysis window ---
win_start <- make_date(1980, 1, 1)
win_end   <- make_date(2018, 12, 1)

# Expect columns: CountryName, SurveyId, PSU, Start.Year, Start.Month, End.Year, End.Month
# Handle missing end dates by treating them as equal to start
exp_prepped <- disas %>%
  mutate(
    Start.Month = coalesce(Start.Month, 1L),
    End.Year    = coalesce(End.Year, Start.Year),
    End.Month   = coalesce(End.Month, Start.Month),
    
    # clamp months into [1..12] just in case
    Start.Month = pmin(pmax(Start.Month, 1L), 12L),
    End.Month   = pmin(pmax(End.Month,   1L), 12L),
    
    start_date0 = make_date(Start.Year, Start.Month, 1),
    end_date0   = make_date(End.Year,   End.Month,   1),
    
    # fix inverted intervals if any (swap)
    start_date  = pmin(start_date0, end_date0),
    end_date    = pmax(start_date0, end_date0),
    
    # clip to the analysis window
    start_clip  = pmax(start_date, win_start),
    end_clip    = pmin(end_date,   win_end)
  ) %>%
  # keep intervals that intersect the window
  filter(start_clip <= end_clip)

# Expand each disaster row to one row per month it spans
exp_months <- exp_prepped %>%
  mutate(month_seq = map2(start_clip, end_clip, ~ seq(.x, .y, by = "month"))) %>%
  dplyr::select(CountryName, SurveyId, psu, month_seq) %>%
  unnest(month_seq) %>%
  mutate(month_label = format(month_seq, "%Y_%m"),
         value = 1L)

# Build the wide month-indicator table (one col per month)
exp_months <- exp_months %>%
  group_by(CountryName, SurveyId, psu, month_label) %>%
  summarise(value = max(value), .groups = "drop") %>%
  pivot_wider(names_from = month_label, values_from = value,
              values_fill = 0)

# Ensure all months from 1980-01 to 2018-12 exist as columns (even if 0)
all_months <- format(seq(win_start, win_end, by = "month"), "%Y_%m")
missing_cols <- setdiff(all_months, names(exp_months))
if (length(missing_cols)) {
  exp_months[missing_cols] <- 0L
}

# Arrange columns in chronological order
exp_months <- exp_months %>%
  dplyr::select(CountryName, SurveyId, psu, all_of(sort(all_months)))

# Convert from wide to long
exp_windows <- exp_months %>%
  pivot_longer(
    cols = matches("^\\d{4}_\\d{2}$"),          # all columns like 1981_01, 1981_02, ...
    names_to = c("year", "month"),              # split into separate columns
    names_sep = "_",                            # split by underscore
    values_to = "disaster"                      # new value column
  ) %>%
  mutate(
    year = as.integer(year),
    month = as.integer(month)
  ) %>%
  arrange(CountryName, SurveyId, psu, year, month)  


### Identify if a disaster occurred in any of the past 3 to 24 months 
library(dplyr)
library(zoo)

# 3-month exposure window
exp_windows <- exp_windows %>%
  arrange(CountryName, SurveyId, psu, year, month) %>%          # chronological order
  group_by(CountryName, SurveyId, psu) %>%
  mutate(
    # Convert year-month to numeric index if needed
    ym_index = (year - min(year)) * 12 + month,
    disaster_3m = rollapply(
      disaster,
      width = 3,             # rolling window length (3 months)
      align = "right",       # rolling over current + previous 2 months
      fill = 0,              # fill first few with 0
      FUN = max,
      partial = TRUE         # allow partial windows at start
    )
  ) %>%
  ungroup()


# 6-month exposure window
exp_windows <- exp_windows %>%
  arrange(CountryName, SurveyId, psu, year, month) %>%          # chronological order
  group_by(CountryName, SurveyId, psu) %>%
  mutate(
    # Convert year-month to numeric index if needed
    ym_index = (year - min(year)) * 12 + month,
    disaster_6m = rollapply(
      disaster,
      width = 6,             # rolling window length (3 months)
      align = "right",       # rolling over current + previous 2 months
      fill = 0,              # fill first few with 0
      FUN = max,
      partial = TRUE         # allow partial windows at start
    )
  ) %>%
  ungroup()


# 9-month exposure window
exp_windows <- exp_windows %>%
  arrange(CountryName, SurveyId, psu, year, month) %>%          # chronological order
  group_by(CountryName, SurveyId, psu) %>%
  mutate(
    # Convert year-month to numeric index if needed
    ym_index = (year - min(year)) * 12 + month,
    disaster_9m = rollapply(
      disaster,
      width = 9,            # rolling window length (3 months)
      align = "right",       # rolling over current + previous 2 months
      fill = 0,              # fill first few with 0
      FUN = max,
      partial = TRUE         # allow partial windows at start
    )
  ) %>%
  ungroup()


# 12-month exposure window
exp_windows <- exp_windows %>%
  arrange(CountryName, SurveyId, psu, year, month) %>%          # chronological order
  group_by(CountryName, SurveyId, psu) %>%
  mutate(
    # Convert year-month to numeric index if needed
    ym_index = (year - min(year)) * 12 + month,
    disaster_12m = rollapply(
      disaster,
      width = 12,            # rolling window length (3 months)
      align = "right",       # rolling over current + previous 2 months
      fill = 0,              # fill first few with 0
      FUN = max,
      partial = TRUE         # allow partial windows at start
    )
  ) %>%
  ungroup()


#### Link with the DHS data

## Recent exposure - past 3 to 12 months from the interview date
data_exp <- data %>%
  left_join(
    exp_windows,
    by = c(
      "CountryName" = "CountryName",
      "SurveyId"    = "SurveyId",
      "psu"         = "psu",
      "intYr"       = "year",
      "intMo"       = "month"
    )
  )

## Early life exposure - infancy, age 1, age 2, in-utero
data_exp <- data_exp %>% 
  mutate(birthYr_plus1 = birthYr+1) %>% 
  mutate(birthYr_plus2 = birthYr+2) %>% 
  mutate(birthYr_plus3 = birthYr+3) 


exp_windows_infancy <- exp_windows %>% 
  dplyr::select(CountryName, SurveyId, psu, year, month, disaster_12m) %>% 
  rename(disaster_infancy = disaster_12m)

exp_windows_age1 <- exp_windows %>% 
  dplyr::select(CountryName, SurveyId, psu, year, month, disaster_12m) %>% 
  rename(disaster_age1 = disaster_12m)

exp_windows_age2 <- exp_windows %>% 
  dplyr::select(CountryName, SurveyId, psu, year, month, disaster_12m) %>% 
  rename(disaster_age2 = disaster_12m)

exp_windows_inutero <- exp_windows %>% 
  dplyr::select(CountryName, SurveyId, psu, year, month, disaster_9m) %>% 
  rename(disaster_inutero = disaster_9m)

data_exp <- data_exp %>%
  left_join(
    exp_windows_infancy,
    by = c(
      "CountryName" = "CountryName",
      "SurveyId"    = "SurveyId",
      "psu"         = "psu",
      "birthYr_plus1" = "year",
      "birthMo"       = "month"))

data_exp <- data_exp %>%
  left_join(
    exp_windows_age1,
    by = c(
      "CountryName" = "CountryName",
      "SurveyId"    = "SurveyId",
      "psu"         = "psu",
      "birthYr_plus2" = "year",
      "birthMo"       = "month"))

data_exp <- data_exp %>%
  left_join(
    exp_windows_age2,
    by = c(
      "CountryName" = "CountryName",
      "SurveyId"    = "SurveyId",
      "psu"         = "psu",
      "birthYr_plus3" = "year",
      "birthMo"       = "month"))

data_exp <- data_exp %>%
  left_join(
    exp_windows_inutero,
    by = c(
      "CountryName" = "CountryName",
      "SurveyId"    = "SurveyId",
      "psu"         = "psu",
      "birthYr" = "year",
      "birthMo"     = "month"))




data_exp <- data_exp %>% 
  # Restrict data to 1980-2018 (period that the geo-referenced disaster data spans)
  filter(intYr>=1980 & intYr<=2018) %>% 
  # Regions where no disasters occurred in the period have NA values - correct to 0
  mutate(disaster = ifelse(is.na(disaster), 0, disaster)) %>% 
  mutate(disaster_3m = ifelse(is.na(disaster_3m), 0, disaster_3m)) %>% 
  mutate(disaster_6m = ifelse(is.na(disaster_6m), 0, disaster_6m)) %>% 
  mutate(disaster_9m = ifelse(is.na(disaster_9m), 0, disaster_9m)) %>% 
  mutate(disaster_12m = ifelse(is.na(disaster_12m), 0, disaster_12m)) %>% 
  mutate(disaster_inutero = ifelse(is.na(disaster_inutero), 0, disaster_inutero)) %>% 
  mutate(disaster_infancy = ifelse(is.na(disaster_infancy), 0, disaster_infancy)) %>% 
  mutate(disaster_age1 = ifelse(is.na(disaster_age1), 0, disaster_age1)) %>% 
  mutate(disaster_age2 = ifelse(is.na(disaster_age2), 0, disaster_age2))
  

saveRDS(data_exp, "Data/data_for_analysis_emdat.rds")

summary <- data_exp %>% 
  group_by(CountryName, disaster_6m) %>% 
  summarise(count = n())