README.Rmd

---
title: "Reproducible materials"
author: "Diego Alburez-Gutierrez"
output:
  github_document:
    toc: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

This document presents the code needed to reproduce the empirical analysis, figures, and tables presented in the paper:

Alburez-Gutierrez, D. (accepted). *The demographic drivers of grief and memory after genocide in Guatemala.* Demography. Temporary DOI: https://dx.doi.org/10.4054/MPIDR-WP-2021-003.

The package `EGM` is needed to reproduce the results of this study. It can be downloaded from Github using the `devtools` package in R:

```{r, eval = FALSE}
  library(devtools)
  install_github("alburezg/EGM", dep = FALSE)
```

After installation, the package can be loaded, together with the other packages needed for the analysis (see version numbers and other technical details at the end of this document): 


```{r}
library(EGM)

library(tidyverse)
library(knitr)
library(kableExtra)
library(ggrepel)
library(patchwork)
```

```{r, eval=T, include=F, warning=F, message=F}
# This private package contains the actual data used to reproduce the figures. 
# This package is used to produce the results in this report, but is not publically available given privacy concerns
library(rionegro)
```

# Structure of the data

The microdata data from Rio Negro cannot be shared in full given privacy concerns. The raw data includes sensitive information such as date of birth, date of death, and kinship ties, which could be used to uniquely identify members of the population. Nevertheless, this section provides a brief overview of the seven datasets referenced in the code below. Note that the datasets included in the EGM package are for illustrative purposes and only include the columns referenced in the analysis.

The `final` dataset records all members of the community, reconstructed using the EGM method (mock data for display purposes only). This is a dataset without duplicates in which individuals are rows. The `ind_id` column is a unique identifier which can be used to locate the relatives of a given individual (e.g., the person with `ind_id == 1` is the daughter of the person with `ind_id == 8`). The column  `id_cuest` indicates in which questionnaires a given individual was reported (see discussion on 'multiple genealogical reporting in the paper'). All other columns are self-explanatory (`DoB` and `DoD` are dates of birth and death; `alive` indicates whether a given person was alive when the data were collected in 2016; `origin_local` is used to identify individuals born in Rio Negro or Pacux).

```{r}
EGM::final %>% select(-DoD_with_alive) %>% head %>% kable
```

# Reproducing results from the paper

In what follows, I show how these datasets were analysed to produce the results reported in the paper. Note that the code will not work using only the sample datasets included in the EGM package. Researchers interested in replicating any particular section of the analysis are encouraged to contact the author. 

## Table 1

### Define parameters

```{r}
war_year <- 1982

death_year_min <- 1983
birth_year_min <- 1982
  
custom_interval <- c(0,15,45)
levels_ordered <- c("0 to 14", "15 to 44", "45+")
levels_custom <-  c("15 to 44", "0 to 14",  "45+")

cos <- war_year - custom_interval

cos_labs <- c(
      paste0(cos[2]+1, "-", cos[1])
      , paste0(cos[3], "-", cos[2])
      , paste0("Pre-",cos[3])
      , "All cohorts"
    )
```

### Define functions

```{r}
format_zero_table <- function(x) gsub("^0\\.", ".", x)

census <- function(df, from, to = NA, subset_df = T, born_in_RN = F, born_bv = F,
                   gen_age = T, women_only = F, men_only = F,
                   by = 5, update_relatives = T, show_messages = F,
                   subset_on_DoB = F, min_DoB, max_DoB) {
  # If subset_df, it returns a DF with only the relevant cases
  # Otherwise, it will create a column 'keep' to indicate which values should be kept.
  if (is.na(to)) to <- from
  # keep only cases where dates are known
  # set DoD of alive people as 2010
  if (from == to) df$census <- from
  else df$census <- paste(from,to,sep = "-")
  df$DoD_temp <- df$DoD
  df$DoD_temp[df$alive == 1] <- 2100

  if (born_in_RN) df <- dplyr::filter(df, origin_literal == 1)

  if (born_bv) df <- dplyr::filter(df, origin_departamento == "baja verapaz")

  if (women_only) df <- df %>% dplyr::filter(sex == 2)

  if (men_only) df <- df %>% dplyr::filter(sex == 1)

  if (subset_df) {
    df <- df[!is.na(df$DoB) & !is.na(df$DoD_temp),]
    DoD <- df$DoD_temp
    df$DoD_temp <- NULL
    DoB <- df$DoB
    # keep only those alive in given interval
    df <- df[DoB <= to & DoD >= from,]
  }
  else {
    DoD <- df$DoD_temp
    df$DoD_temp <- NULL
    DoB <- df$DoB
    df$keep <- NA
    df$keep <- ifelse(DoB <= to & DoD >= from, T, F)
  }
  if(subset_on_DoB) {
    # keep only those alive in given interval
    df <- df %>%
      dplyr::filter(DoB >= min_DoB & DoB <= max_DoB)
  }
  if (gen_age) {
    df$age_at_census <- round(mean(c(to,from))) - df$DoB
    df$age_gr_at_census <- age_group_custom(df$age_at_census, by = by, as_factor = T)
  }
  if (update_relatives) {
    # keep original copy
    df$children_original <- df$children
    df$father_original <- df$father
    df$mother_original <- df$mother
    df$parents_original <- df$parents
    df$cousins_original <- df$cousins
    df$grandparents_original <- df$grandparents
    df$granchildren_original <- df$grandchildren
    df$spouses_original <- df$spouses
    if (show_messages) {
      print("Columns substituted for all relative types.")
      print("Original values were saved in '[relative]_original' columns.")
      print("For spouses: Alive during census and older than 14.")
      print("For everyone else, alive during census.")
    }
  }
  return(df)
}

```

### Run analysis

```{r}

 c_df <- 
    final %>%
    dplyr::filter(died_conflict & DoB <= birth_year_min) %>%
    dplyr::mutate(
      sex = sex2,
      age_conflict = ifelse(DoB <= war_year, war_year-DoB, NA),
      ag_conflict = cut(age_conflict, c(custom_interval, Inf), right = F, labels = F),
      ag_conflict = factor(plyr::mapvalues(ag_conflict, from = seq_along(levels_ordered), 
                                           to = levels_ordered), levels = levels_ordered)
    ) %>%
    group_by(ag_conflict, sex) %>%
    dplyr::summarise(n = n())
  
# For census data
  
  all_df <-
    census(final, from = 1980, to = 1984) %>% 
    dplyr::filter(DoB <= birth_year_min) %>%
    dplyr::mutate(
      sex = sex2,
      age_conflict = ifelse(DoB <= war_year, war_year-DoB, NA),
      ag_conflict = cut(age_conflict, c(custom_interval, Inf), right = F, labels = F),
      ag_conflict = factor(plyr::mapvalues(ag_conflict, from = seq_along(levels_ordered), 
                                           to = levels_ordered), levels = levels_ordered)
    ) %>%
    group_by(ag_conflict, sex) %>%
    dplyr::summarise(n = n())
  
  df <- merge(c_df, all_df, by = c("ag_conflict", "sex"))
  
  
  df$share <- (as.numeric(df$n.x)/df$n.y)
  
  df$ag_conflict <- factor(df$ag_conflict, levels = levels(all_df$ag_conflict))
  
  df2 <-
    df %>%
    select(ag_conflict, sex, b.total_original = n.y, c.total_died = n.x, d.share = share) %>%
    reshape(., idvar = "ag_conflict", timevar = "sex", direction = "wide") %>%
    select(
      ag_conflict, 
      original_female = b.total_original.Female, 
      original_male = b.total_original.Male,
      died_female = c.total_died.Female,
      died_male = c.total_died.Male,
      share_female = d.share.Female,
      share_male = d.share.Male
    ) %>%
    dplyr::mutate(
      original_all = original_female + original_male,
      died_all = died_female + died_male,
      share_all = rowMeans(.[ , c("share_female", "share_male") ] )
    ) %>%
    select(ag_conflict, original_all , original_female , original_male,
           died_all , died_female , died_male,
           share_all , share_female , share_male) %>%
    mutate_if(is.numeric, round, 2) %>%
    arrange(ag_conflict) %>%
    dplyr::mutate(ag_conflict = as.character(ag_conflict)) %>%
    data.frame
  
  newrow <- c("Total",colSums(df2[,2:ncol(df2)]))
  
  df3 <- rbind( df2,
                rbind(df2, newrow) %>%
                  dplyr::filter(ag_conflict == "Total") %>%
                  mutate_if(is.character, as.numeric) %>%
                  dplyr::mutate(
                    ag_conflict = "Total",
                    share_all = died_all / original_all *1,
                    share_female = died_female / original_female*1,
                    share_male = died_male / original_male*1
                  ) %>%
                  mutate_if(is.numeric, round, 2) %>%
                  data.frame
  )

  em_df <- 
  data.frame(cos_labs) %>% 
  bind_cols(df3) %>% 
  mutate(ag_conflict = dplyr::recode(ag_conflict, "Total" = "All ages")) %>% 
  mutate(
    across(where(is.numeric), round, 2)
    , across(where(is.numeric), as.character)
    , across(everything(), format_zero_table)
  )
  
  em_df %>% 
    kable(
      format = "simple"
      , booktabs = T, label = "excess_mort"
      , col.names = c("","", rep(c("All", "Female", "Male"), 3))
      , caption = "Age-sex structure of the Rio Negro population before the 1982 killings and excess deaths from the killings. (A) Population by age and sex at the start of 1982; (B) direct deaths from the killings; (C) proportion of the pre-genocide population killed."
    ) %>% 
    kable_styling() %>% 
    add_header_above(c(
      "", "", "(A)" = 3
      , "(B)" = 3
      , "(C) = B/A" = 3
    )) %>% 
    add_header_above(c(
      "Cohort" = 1
      , "Age in 1982" = 1
      , "Alive in 1982" = 3
      , "Killed (total)" = 3
      , "Proportion killed" = 3
    ), line = F) 
  
write.csv(em_df, "Tab1.csv", row.names = F)
```

## Table 2

### Define parameters

```{r}

br <- c(0, 15, 45, Inf)

rel_old <- c(
  "children", "parents", "siblings", "spouses"
  , "grandchildren", "grandparents", "cousins"
  , "nuclear", "extended"
  )


rel_new <- c(
  "Children", "Parents", "Siblings", "Spouses"
  , "Grandchildren", "Grandparents", "Cousins"
  , "Nuclear (any)", "Extended (Total)"
) 

sex_ag_levs <- 
  c("Female_[0,15]", "Female_(15,45]", "Female_(45,Inf]", "Female_Total"
    , "Male_[0,15]", "Male_(15,45]", "Male_(45,Inf]", "Male_Total"
    , "Total_Total"
    )


```

### Define functions

```{r}
create_age_groups <- function(df, col_in, br = c(0, 15, 45, Inf), include.lowest = T, col_out = "ag"){
  new <- cut(df[, col_in], br, include.lowest = include.lowest)
  df[, col_out] <- new
  df
}
```

### Run analysis

```{r}

expo_all <- 
  census(final, 1983) %>% 
  create_age_groups(col_in = "age_at_census", br = br) %>% 
  select(
    sex = sex2, ag, nuclear = c_index_nuclear, extended = c_index_extended
    , ends_with("_death"), -age_death
  ) %>% 
  pivot_longer(-c(sex, ag)) %>% 
  mutate(name = gsub("_death", "", name))

# Fix error
expo_all$value[expo_all$ag == "[0,15]" & expo_all$name == "grandchildren"] <- 0

# SUmmary by age group
expo <-
  expo_all %>% 
  group_by(sex, ag, name) %>% 
  summarise(
    sd = sd(value, na.rm = T)
    , value = mean(value, na.rm = T)
    ) %>% 
  ungroup() %>% 
  mutate(
      pretty = paste0(round(value, 1), " (", round(sd, 1), ")")
    ) %>% 
  select(-sd)

# Summary for all ages
expo_all_ag <-
  expo_all %>% 
  group_by(sex, name) %>% 
  summarise(
    sd = sd(value, na.rm = T)
    , value = mean(value, na.rm = T)
    ) %>% 
  ungroup() %>% 
  mutate(
    # pretty = paste0(hide0(round(value, 1)), " (", hide0(round(sd, 1)), ")")
        pretty = paste0(round(value, 1), " (", round(sd, 1), ")")
    , ag = "Total"
    ) %>% 
  select(-sd)

# Summary for all ages and sexes
expo_all_ag_sex <-
  expo_all %>% 
  group_by(name) %>% 
  summarise(
    sd = sd(value, na.rm = T)
    , value = mean(value, na.rm = T)
    ) %>% 
  ungroup() %>% 
  mutate(
    # pretty = paste0(hide0(round(value, 1)), " (", hide0(round(sd, 1)), ")")
        pretty = paste0(round(value, 1), " (", round(sd, 1), ")")
      , sex = "Total"
    , ag = "Total"
    ) %>% 
  select(-sd)

expo_df <- 
  expo %>% 
  bind_rows(expo_all_ag) %>% 
  bind_rows(expo_all_ag_sex) %>% 
  mutate(
    name = plyr::mapvalues(name, rel_old, rel_new)
    , name = factor(name, levels = rel_new)
  ) %>% 
  select(-value) %>%
  pivot_wider(names_from = c(sex, ag), values_from = pretty) %>% 
  arrange(name)  

expo_df2 <- expo_df[ , c("name", sex_ag_levs)] 

expo_df2 %>% 
  kable(
    format = "pipe"
      , booktabs = T
    , caption = "Number of relatives lost by sex and age of genocide survivors at the time of the killings (mean and standard deviation)."
    , label = "kin_loss_mean"
    , col.names = c(" ", "F 0-14", "F 15-44", "F 45+", "F Tot",  "M 0-14", "M 15-44", "M 45+", "M Tot", "All_Total")
  ) 

write.csv(expo_df2, "Tab2.csv", row.names = F)

```

## Table 3

### Define parameters

```{r}
deaths_keep <- 1:5

col_kin <- c(
  "nuclear (any)"
  , "children", "parents", "siblings", "spouses"
  , "extended (any)"
  , "grandchildren","grandparents", "cousins"
    )

```


## Run analysis

```{r}
pop_temp <- 
  census(final, 1982) %>%
  dplyr::filter(!died_conflict) %>%
  mutate(
    died_other = ifelse(alive == 0 & !died_conflict, T, F)
    # THis creates some NAs, presumably for people whose current state is unknown
    # Assume they are dead, for convenience
    , died_other = ifelse(is.na(died_other), T, F)
  )


prog_df <-
  pop_temp %>% 
  select(
    tidyselect::ends_with("_death")
    , `extended (any)` = c_index_extended, `nuclear (any)` = c_index_nuclear
    ) %>% 
  select(contains(col_kin)) %>% 
  pivot_longer(everything(), names_to = "kin", values_to = "deaths") %>% 
  mutate(
    kin = gsub("_death", "", kin)
    # Some issue with codinf for parents
    , deaths = ifelse(is.na(deaths), 0, deaths)
    ) %>% 
  count(kin, deaths) %>% 
  tidyr::complete(kin, deaths, fill = list(n = 0)) %>% 
  filter(!is.na(deaths)) %>% 
  group_by(kin) %>% 
  # n is number of poeople who had at least x deaths
  mutate(
    # People who experienced x or more deaths
    # This is estimated like Tx in a life table
    risk = rev(cumsum(rev(n))) 
    , prog = risk / lag(risk)
    , prog = ifelse(is.nan(prog), NA, prog)
    , percent = round(prog * 100, 1)
  ) %>% 
  ungroup() %>% 
  data.frame()

prog_temp <- 
  prog_df %>% 
  filter(deaths %in% deaths_keep) %>%
  select(kin, deaths, risk, prob = prog) %>% 
  mutate(
    kin = factor(kin, levels = col_kin)
    , deaths = paste0(deaths-1, " to ", deaths, "+")
    , prob = ifelse(kin == "parents" & deaths == 3, NA, prob)
    , prob = ifelse(is.na(prob), "", round(prob, 2))
    , out = paste0(prob, " [", risk, "]")
    , out = gsub("^0\\.", ".", out)
    ) %>% 
  filter(prob > 0) %>% 
  arrange(kin) %>% 
  pivot_wider(-c(risk, prob), names_from = deaths, values_from = out, values_fill = list(out = "")) %>% 
  mutate(kin = str_to_title(kin)) %>% 
  rename(`(n-1) to n+` = kin) 

prog_temp %>% 
  kable(format = "simple"
    , booktabs = T, escape = F
    , caption = "Proportion of genocide survivors who experienced the death of $n$ or more relatives ($n+$), conditional on having experienced the death of at least $(n-1)$ relatives. In square brackets, population who experienced at least $n$ deaths of a given relative type."
    , label = "bpr"
    ) 

write.csv(prog_temp, "Tab3.csv", row.names = F)

```

## Figure 1

### Define parameters

```{r}

# For data wrangling

awareness_ages <- c(5, 12)

# For plot

y_br_small <- c(1982, 1990, 2000, 2010, 2015)

min_br <- seq(1980, 2015, 5)

l1 <- "Genocide Witnesses\n(proportion of population)"
l2 <- "Mean Age"

lab_df <- data.frame(
  year = c(1995, 2005)
  , value = c(14, 40)
  , type = c(l1, l2)
)

```

### Define functions

```{r}
# Given a type of relative and awareness age, determine the share
# of witnesses in a populatio over time, or the share of pop with 
# a victim from the genocide
get_witness_df <- function(family_type, awareness_age){
  
  # Parameters
  relatives <- switch (family_type,
                       nuclear = c("father", "mother", "siblings", "children", "spouses")
                       , extended = c("father", "mother", "siblings", "children", "spouses", "cousins", "grandparents", "grandchildren", "uncles")
                       , parents = c("father", "mother")
  )
  
  year_range <- 1982:2016
  
  death_year_min <- 1983
  birth_year_min <- 1982
  
  war_year <- 1982
  
  age_br <- c(0, 15, 45, 100)
  
  birth_year_witness_min <- birth_year_min - awareness_age
  
  # 1. Get data
  
  everyone <- 
    final %>%
    select(
      ind_id, DoB, DoD, DoD_with_alive,alive, died_conflict, sex = sex2
      , father, mother, siblings, children, spouses
      , cousins, grandparents, grandchildren, uncles
    ) %>% 
    dplyr::mutate(
      age_conflict = ifelse(DoB <= war_year, war_year-DoB, NA),
      ag_conflict = cut(age_conflict, age_br, include.lowest = T)
    )
  
  # 1.1 First-hand witnesses
  
  alive_1982 <- 
    everyone %>% 
    dplyr::filter(DoB <= war_year) %>%
    dplyr::filter(DoD_with_alive >= war_year)
  
  # First-hand witness: alive in 1982 and older than awareness_age
  # years old at the time
  
  # Men and women
  witness1_ids <-
    alive_1982 %>% 
    dplyr::filter(DoB < birth_year_witness_min) %>%
    pull(ind_id)
  
  # Assign col to df to idenfify witnesses
  everyone$witness1 <-  everyone$ind_id %in% witness1_ids
  
  # Since weare ignoreing witness 2 here
  everyone$witness2 <-  F
  
  # Need a df BY YEAR
  witness_share <- get_pop_composition_witness(everyone,year_range)
  
  levs <- c("w1_w", "w2_w", "other_w")
  labs_wit <- c("No", "Second-hand", "First-hand")
  y_br <- c(1982, seq(1985, 2015, 5))
  
  share_df <- 
    witness_share %>% 
    pivot_longer(-year) %>% 
    filter(name %in% levs) %>% 
    mutate(name = factor(name, levels = rev(levs)))
  
  # 3 Mean age of witnesses
  
  # 3.1. Men and women together
  
  # Need a df BY YEAR
  witness_age <- get_mean_age_witness(everyone,year_range)
  
  levs_age <- c("w1", "w2", "other")
  labs_wit <- c("No", "Second-hand", "First-hand")
  y_br <- c(1982, seq(1985, 2015, 5))
  
  age_df <- 
    witness_age %>% 
    pivot_longer(-year) %>% 
    mutate(
      status = gsub("_mean|_sd", "", name)
      , measure = gsub("w1_|w2_|other_", "", name)
    )  %>% 
    select(-name) %>% 
    pivot_wider(names_from = measure, values_from = value) %>% 
    mutate(
      low = mean - sd*0.5
      , high = mean + sd*0.5
    ) %>% 
    mutate(status = factor(status, levels = rev(levs_age))) 
  
  both_survivors <- 
    bind_rows(
      share_df %>%
        mutate(
          status = gsub("_w", "", name)
          , status = factor(status, levels = rev(levs_age))
          , type = "Share of population"
        ) %>% 
        select(year, status, value, type)
      , age_df %>% 
        select(year, status, value = mean, age_low = low, age_high = high) %>% 
        mutate(
          type = "Mean Age"
          , status = factor(status, levels = rev(levs_age))
        )
    ) %>% 
    mutate(type = factor(type, levels = c("Share of population", "Mean Age")))
  
  return(both_survivors)
}

# For plots of witness survival
get_pop_composition_witness <- function(everyone, year_range){
  
  out <- 
    lapply(year_range, function(y, everyone){
      df <- census(everyone, y)
      # unweighted
      w1 <- sum(df$witness1) 
      w2 <- sum(df$witness2) 
      other <- sum(!df$witness1 & !df$witness2) 
      # weighted
      w1_w <- w1 / nrow(df)
      w2_w <- w2 / nrow(df)
      other_w <- other / nrow(df)
      # return
      data.frame(year = y, w1, w2, other, w1_w, w2_w, other_w)
    }, everyone)  %>% 
    bind_rows()
  out  
}

# For plots of witness survival
get_mean_age_witness <- function(everyone, year_range){
  
  df_temp <-
    everyone %>% 
    mutate(
      status = ifelse(witness1, "w1", "other")
      , status = ifelse(witness2, "w2", status)
    )
  
  out <- 
    lapply(year_range, function(y, df_temp){
      
      census(df_temp, y) %>% 
        group_by(status) %>% 
        summarise(
          mean = mean(age_at_census, na.rm = T)
          , sd = sd(age_at_census, na.rm = T)
        ) %>% 
        ungroup() %>% 
        pivot_longer(-status) %>% 
        mutate(
          id = paste0(status,"_",name)
          , value = ifelse(is.na(value), 0, value)
        ) %>% 
        select(id, value) %>% 
        pivot_wider(names_from = id) %>% 
        mutate(year = y)
    }, df_temp)  %>% 
    bind_rows()
  out  
}
```

### Run analysis

```{r}

both_survivors <- 
  get_witness_df(
  # Relative type to consider for estimating relatedness to victim
  family_type = "nuclear"
  # Minimum age at which we assume that survivors of the killings still 
  # remember
  , awareness_age = awareness_ages[2]
  ) 

```

### Faceted plot

```{r}

# Share
p_s <- 
  both_survivors %>% 
  filter(
    type == "Share of population"
    , status %in% c("w1")
  ) %>% 
  mutate(
    type = recode(
      type,  "Share of population" = l1, "Mean Age" = l2
    )
    , value = ifelse(type == l1, value*100, value)
    , value = round(value/100, 2)
  ) %>% 
    filter(between(year, 1983, 2015)) %>% 
  ggplot() +
  # Age
  geom_line(
    aes(x = year, y = value, linetype = status)
    , show.legend = F
    , size = 2
  ) +
  geom_vline(xintercept = 1982, linetype = "dashed") +
  scale_y_continuous(l1) +
  scale_x_continuous("Year", breaks = y_br_small, minor_breaks = min_br) +
  scale_color_manual("", values = c("black", "red"), guide=FALSE) + 
  scale_fill_manual("", values = c("black", "red"), guide=FALSE) + 
    scale_linetype_manual("", values = c(1:2), labels = c("Witness", "Nonwitness")) +
  coord_cartesian(ylim = c(0.05, 0.65)) + 
  theme_bw(base_size = 18) 

# Age
p_a <- 
  both_survivors %>% 
  filter(type == "Mean Age") %>% 
  mutate(
    type = recode(
      type,  "Share of population" = l1, "Mean Age" = l2
    )
    , value = ifelse(type == l1, value*100, value)
    , status = factor(status, levels = c("w1", "w2", "other"))
  ) %>% 
  filter(between(year, 1983, 2015)) %>% 
  ggplot() +
  # Age
  geom_line(
    aes(x = year, y = value, linetype = status)
    # , show.legend = F
    , size = 2
  ) +
  # SurvWitness age bands
  geom_ribbon(
    aes(x = year, y = value, ymin = age_low, ymax = age_high)
    , alpha = 0.3, show.legend = F
    , data = . %>% filter(status %in% c("w1"))
  ) +
  # Nonwitness age bands
  geom_ribbon(
    aes(x = year, y = value, ymin = age_low, ymax = age_high)
    , alpha = 0.3, show.legend = F
    , data = . %>% filter(status %in% c("other"))
  ) +
  geom_vline(xintercept = 1982, linetype = "dashed") +
  scale_y_continuous(l2) +
  scale_x_continuous("Year", breaks = y_br_small, minor_breaks = min_br) +
  scale_color_manual("", values = c("black", "red"), guide=FALSE) + 
  scale_fill_manual("", values = c("black", "red"), guide=FALSE) + 
  # scale_linetype_manual("", values = c(2:1), labels = c("Nonwitness","Witness")) +
  scale_linetype_manual("", values = c(1:2), labels = c("Witness", "Nonwitness")) +
  coord_cartesian(ylim = c(5, 65)) + 
  theme_bw(base_size = 18) +
  guides(linetype = guide_legend(override.aes = list(size = 1))) 

combined <- 
  p_s + p_a & theme(
  legend.position = "bottom"
  , legend.key.width = unit(2,"cm")
  , axis.title.y = element_text(face="bold")    
  , axis.title.x = element_text(face="bold")
  )

combined + plot_layout(guides = "collect")

ggsave(paste0("Fig1.pdf"), width = 25, height = 12, units = "cm")

```

## Figure 2

What does 'overlapped with victim' mean? 

```{r}



# rcode ldh129
data <- 
  data.frame(
    y1 = 1:5
    , y2 = 1:5
    , x1 = c(1975, 1975, 1980, 2010, 1990)
    , x2 = c(1982, 1982, 2000, 2019, 2019)
    , shape_start = c("19","19","19", "19", "19")
    , shape_end = c("4", "4", "4",NA, NA)
    , colour = c("red", "red", "black", "black", "black")
  ) 
  
data %>% 
  ggplot()  +
  geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2, colour = colour), size = 1) +
  # start points
  geom_point(aes(x = x1, y = y1, shape = shape_start, colour = colour), size = 4) +
  # end points
  geom_point(aes(x = x2, y = y1, shape = shape_end, colour = colour), size = 4) + 
  geom_vline(xintercept = 1982, linetype = "dashed") + 
  scale_colour_manual(values = c("black", "red"), guide = F) +
  scale_shape_manual(values = c(19, 4), guide = F) + 
  # labels
  geom_label_repel(
    aes(x = x, y = y, label = label)
    , nudge_x = 10
    , nudge_y = 0.8
    , segment.curvature = 0.1
    , segment.ncp = 3
    , segment.angle = 20
    , arrow = arrow(length = unit(0.05, "npc"))
    , data = data.frame(
      x = c(1981)
      , y = c(3.1)
      , label = c("Overlapped with victim")
    )
  ) +
  scale_y_continuous("Sibling ID", breaks = 1:5, labels = LETTERS[1:5], minor_breaks = 1:5) +
  scale_x_continuous("Year", breaks = c(1982, 1990, 2000, 2010, 2015), minor_breaks = min_br) +
  coord_cartesian(xlim = c(1980,2015), ylim = c(0.5, 5.5)) +
  theme_bw() +
  theme(
   axis.title.y = element_text(face="bold")    
  , axis.title.x = element_text(face="bold")
  )

ggsave("Fig2.pdf", width = 15, height = 5, units = "cm")

```

## Figure 3

### Define parameters

```{r}
war_year <- 1982

age_br <- c(0, 15, 45, 100)
y_br_small <- c(1980, 1990, 2000, 2010, 2015)

nuc <- c("parents", "siblings", "children", "spouses")
fam_levels <- c("Extended family", "Nuclear family")
knew_br <- c("met_relative","all", "did_not_meet_relative")
knew_levs <- c("Overlapped with a victim","Related to a victim", "No overlap")

```

### Define functions

```{r}
get_df_share_bereaved_over_time <- function(condition = "all"){
  
  # General dfs
  
  everyone <- 
    final %>%
    select(
      ind_id, DoB, DoD, DoD_with_alive,alive, died_conflict, sex = sex2
      , father, mother, siblings, children, spouses
      , cousins, grandparents, grandchildren, uncles
    ) %>% 
    dplyr::mutate(
      age_conflict = ifelse(DoB <= war_year, war_year-DoB, NA),
      ag_conflict = cut(age_conflict, age_br, include.lowest = T)
    )
  

  # In order to determine whether the populations actually knew
  # (ie coexisted) with dead realtives, some filtering is required.
  # The trick is that these measure show all refer to the same
  # population, which is why filtering on population by date of
  # birth in here is unwise.
  year_range <- 1982:2016
  
  baseline_pop <- 
    everyone %>% 
    dplyr::filter(DoD_with_alive >= war_year) %>% 
    # add columns for each condition
    mutate(
      all = T
      , met_relative = ifelse(DoB <= war_year & DoD_with_alive >= war_year, T, F)
      # , did_not_meet_relative = ifelse(DoB > 1983 & DoD_with_alive >= war_year, T, F)
      , did_not_meet_relative = !met_relative
    )
  
  # Get all relative types
  rels <- list(
    parents = c("father", "mother"),
    siblings = "siblings", 
    children = "children", 
    spouses = "spouses", 
    cousins = "cousins", 
    grandparents = "grandparents", 
    grandchildren = "grandchildren", 
    uncles = "uncles",
    nuclear = c("father", "mother", "siblings", "children", "spouses"),
    extended = c("father", "mother", "siblings", "children", "spouses", "cousins", "grandparents", "grandchildren", "uncles")
  )
  
  
  share_df <- 
    lapply(seq_along(rels), function(n,condition){
      # print(n)
      col_name <- names(rels)[n]
      
      baseline_pop$has_a_relative_victim <- 
        check_if_any_relative_fulfils_condition(
          df = baseline_pop
          , ref = everyone
          , status_col = "died_conflict"
          , relatives = rels[[n]]
        )
      
      pop_share <- 
        get_pop_composition_relative_died(
          everyone = baseline_pop
          , year_range
          , col_name = col_name
          , condition
          )
      
      pop_share
      
    }, condition) %>% 
    bind_rows()
  # reduce(left_join, "year")
  
}

# FOr witness analysis
check_if_any_relative_fulfils_condition <- function(df, ref, status_col, relatives = c("father", "mother", "siblings", "children", "spouses")){
  # Create a column with all relatives pasted
  df_witnesses <- 
    df %>% 
    mutate(relatives = find_relatives(. , relatives)) %>% 
    select(ind_id, relatives, status = starts_with(status_col))
  
  # Now, loop over each element and find people who match the given 
  # status col criteria
  any_relative_fulfils_condition <- 
    lapply(df_witnesses$relatives, function(rel, ref){
      ids <- unlist(strsplit(rel, ";"))
      ids <- as.numeric(ids[!(ids == "NA"|is.na(ids))])
      relatives_status <- ref[ , status_col][match(ids,ref$ind_id )]  
      any(relatives_status)
    }, ref) %>% 
    unlist()
  return(any_relative_fulfils_condition)
}

# For any RN df, returns a vector with all requested relatives
find_relatives <- function(df, relatives){
  # Create a vector with all relevant ids of relative
  rel_df <- df[ , relatives] 
  if(!is.data.frame(rel_df)) rel_df <- as.data.frame(rel_df) # added on 15.10.2020
  rels <- NULL
  
  for(row in 1:nrow(rel_df)) {
    new <- unlist(rel_df[row,])
    # new[new == "NA"] <- NA
    new <- new[!is.na(new)]
    new_vec <- paste0(new, collapse = ";")
    rels[row] <- new_vec
  }
  
  rels[rels == ""] <- NA
  return(rels)
}

get_pop_composition_relative_died <- function(everyone, year_range, col_name, condition){
  
  out <- 
    lapply(year_range, function(y, everyone, condition){
      df <- census(everyone, y)
      # total number of bereaved people
      # filtered by condition!!
      pop <- sum(df[df[,condition], "has_a_relative_victim"]) 
      # weighted
      share <- pop / nrow(df)
      # return
      data.frame(year = y, pop, share)
    }, everyone, condition)  %>% 
    bind_rows() %>% 
    mutate(relative = col_name)
  out  
}


```

### Run analysis

```{r}
share_are_bereaved_met <- 
  get_df_share_bereaved_over_time(condition = knew_br[1])

share_are_bereaved_all <- 
  get_df_share_bereaved_over_time(condition = knew_br[2])

share_are_bereaved_not_met <- 
  get_df_share_bereaved_over_time(condition = knew_br[3])

share_are_bereaved <- bind_rows(
   share_are_bereaved_met %>% mutate(knew = knew_levs[1])
   , share_are_bereaved_all %>% mutate(knew = knew_levs[2])
  , share_are_bereaved_not_met %>% mutate(knew = knew_levs[3])
) %>% 
  mutate(
    type = ifelse(relative %in% c(nuc, "nuclear"), "Nuclear family", "Extended family")
    , relative = recode(
      relative, 
      "nuclear"="nuclear (any)", "extended"="extended (any)"
      , "uncles" = "aunts/uncles"
      )
    , knew = factor(knew, levels = rev(knew_levs))
    , type = factor(type, levels = fam_levels)
  ) 

```

### Plot

```{r}

share_are_bereaved %>% 
  filter(between(year, 1983, 2015)) %>% 
  filter(knew != knew_levs[3]) %>%
  ggplot() + 
  geom_line(
    aes(x = year, y = share, colour = relative)
    , size = 1
    , show.legend = F
    ) +
  geom_label_repel(
    aes(x = year, y = share, label = relative, colour = relative)
    , show.legend = F
    , data = . %>% filter(year == 2000, knew == knew_levs[2])
  ) +
  geom_vline(xintercept = 1982, linetype = "dashed") +
  scale_y_continuous("Proportion of Population") +
  scale_x_continuous("Year", breaks = y_br_small, labels = y_br_small, minor_breaks = min_br) +
  facet_grid(factor(knew, levels = rev(knew_levs[-3])) ~ type, scales = "fixed") +
  coord_cartesian(xlim = c(1980, 2015)) +
  theme_bw(base_size = 16) +
  theme(
   axis.title.y = element_text(face="bold")    
  , axis.title.x = element_text(face="bold")
  , axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)
  , strip.text.x = element_text(face="bold")
  , strip.text.y = element_text(face="bold")
  )

ggsave("Fig3.pdf", width = 18, height = 15, units = "cm")
```

## Figure 4

### Define parameters

```{r}
wit_br <- c("witness_met", "witness_not_met", "not_witness")
wit_labs <- c("Related and overlapped", "Related but never overlapped", "Not related")

min_age <- 0
max_age <- 75

br_y <- c(seq(0, max_age, 20), max_age)
labs_y <- br_y
labs_y[length(labs_y)] <- paste0(labs_y[length(labs_y)], "+")

filter_years <- c(1983, 2015)
```

### Define functions:

```{r}
# Define function to get ggplot pyramid
get_witness_pyramid <- function(filter_year, wit_br, wit_labs, min_age, max_age){
  
  cell_scale <- 1
  war_year <- 1982
  
  # General dfs
  
  everyone <- 
    final %>%
    select(
      ind_id, DoB, DoD, DoD_with_alive,alive, died_conflict, sex = sex2
      , father, mother, siblings, children, spouses
      , cousins, grandparents, grandchildren, uncles
    ) 
  
  
  # In order to determine whether the populations actually knew
  # (ie coexisted) with dead realtives, some filtering is required.
  # The trick is that these measure show all refer to the same
  # population, which is why filtering on population by date of
  # birth in here is unwise.
  year_range <- 1982:2015
  
  baseline_pop <- 
    everyone %>% 
    dplyr::filter(DoD_with_alive >= war_year) %>% 
    # add columns for each condition
    mutate(
      all = T
      , met_relative = ifelse(DoB <= war_year & DoD_with_alive >= war_year, T, F)
      , did_not_meet_relative = !met_relative
    )
  
  # Get all relative types
  rels <- list(
    parents = c("father", "mother"),
    siblings = "siblings", 
    children = "children", 
    spouses = "spouses", 
    cousins = "cousins", 
    grandparents = "grandparents", 
    grandchildren = "grandchildren", 
    uncles = "uncles",
    nuclear = c("father", "mother", "siblings", "children", "spouses"),
    extended = c("father", "mother", "siblings", "children", "spouses", "cousins", "grandparents", "grandchildren", "uncles")
  )
  
  baseline_pop$has_a_relative_victim <- 
    check_if_any_relative_fulfils_condition(
      df = baseline_pop
      , ref = everyone
      , status_col = "died_conflict"
      , relatives = rels[['extended']]
    )
  
  vic_rel <- 
    baseline_pop %>% 
    select(ind_id, DoB, DoD, DoD_with_alive, sex, has_a_relative_victim, met_relative, did_not_meet_relative) %>% 
    mutate(sex = ifelse(sex == "Female", 2, 1))
  
  witness_df <- 
    vic_rel %>% 
    mutate(
      type = ifelse(!has_a_relative_victim, wit_br[3], NA)
      , type = ifelse(has_a_relative_victim & met_relative, wit_br[1], type)
      , type = ifelse(has_a_relative_victim & did_not_meet_relative, wit_br[2], type)
      , type = factor(type, levels = wit_br)
    )
  
  # To Plot 
  
  pyramid_datapoly_witness_1x1(
    df = witness_df
    , filter_year = filter_year
    , cell_scale = 1
    , min_age
    , max_age
  ) %>% 
    mutate(year = filter_year)
  
}

pyramid_datapoly_witness_1x1 <- function(df, filter_year, cell_scale, min_age, max_age){
  
  by <- cell_scale
  sep <- "-"
  
  cuts2 <- c(seq(min_age,max_age,by))
  
  ag_names <- c( paste0(min_age, sep, by-1), 
                 paste( seq(by,max_age - by, by), seq(by + by - 1,max_age, by), sep =  sep), 
                 paste0(max_age, "+"))
  
  if(by == 1) ag_names <- cuts2
  
  # df of pre-killing population
  pre <- 
    census(final, filter_year) %>% 
    mutate(
      years_survived = DoD_with_alive - filter_year
      , female = ifelse(sex == 1, F, T)
      , age_labs = cut(age_at_census, c(cuts2, Inf), right = F, labels = ag_names)
      , age  = cut(age_at_census, c(cuts2, Inf), right = F, labels = F)
    ) %>% 
    select(id = ind_id, female, DoB, age_labs, age) %>% 
    left_join(df, by = c("id" = "ind_id")) %>% 
    rename(measure = type) %>% 
    mutate(id = as.character(id))
  
  pre_coord <- 
    pre %>% 
    group_by(female, y = age) %>% 
    arrange(measure) %>% 
    mutate(x = 1:n()) %>% 
    ungroup() %>% 
    arrange(y) %>% 
    select(id, female, x, y, measure)
  
  # Get coordinates
  coords <- data.frame(do.call(rbind, apply(pre_coord, 1, pyramid_coords) ))
  
  datapoly <- 
    left_join(pre, coords, by = "id") %>% 
    mutate(
      id = as.character(id)
      # Make women negative
      , x = ifelse(female, -x, x) 
    ) %>% 
    select(-age) %>% 
    rename(age = age_labs)
  
  return(datapoly)
  
}

pyramid_coords <- function(row){
  
    id <- row['id']
    n <- as.numeric(row[c("x", "y")])
    names(n) <- c("x", "y")
    
    x <- c(n['x'], n['x'] + 1, n['x'] + 1, n['x'])
    y <- c(n['y'], n['y'], n['y'] + 1, n['y'] + 1)
    
    data.frame(
      id = id
      , x
      , y
      , row.names = NULL, stringsAsFactors = F)
    
}

```

### Run analysis

```{r}

datapoly <- 
  get_witness_pyramid(filter_years[1], wit_br, wit_labs, min_age, max_age) %>% 
  bind_rows(
    get_witness_pyramid(filter_years[2], wit_br, wit_labs, min_age, max_age)
  )

sex_df <- 
  data.frame(
    x = c( min(datapoly$x ) + sd(datapoly$x), max(datapoly$x ) - sd(datapoly$x) )
    , y = max(datapoly$y - 5)
    , label = c("Female", "Male")
  )

```

### Plot

```{r}
datapoly %>%
  mutate(y = y -1) %>% 
  ggplot(aes(x=x, y=y)) + 
  geom_polygon(
    aes(fill=measure, group=id)
    , colour = "black", size = 0.05
  ) +
  geom_label(aes(x = x, y = y, label = label), data = sex_df) +
  scale_x_continuous("Population"
                     , breaks = scales::pretty_breaks(n = 6)
                     , labels = function(br) abs(br)
  ) +
  scale_y_continuous("Age"
                     , breaks = br_y
                     , labels = labs_y
                     , sec.axis = dup_axis(name = "")
  ) +
  scale_fill_discrete("", breaks = wit_br, labels = wit_labs) +
  facet_grid(~year) +
  coord_fixed() +
  theme_bw() +
    theme(
   axis.title.y = element_text(face="bold")    
  , axis.title.x = element_text(face="bold")
  , strip.text.x = element_text(face="bold")
  , legend.position = "bottom"
  )

ggsave("Fig4.pdf", width = 15, height = 10, units = "cm")
```

# Session information

Report by Diego Alburez-Gutierrez - https://twitter.com/d_alburez

```{r, echo=FALSE, warning=FALSE}
print(paste("Report created:", Sys.time()))
sessionInfo()
```