alburezg/genocide_guatemala

Reproducible materials or paper: Alburez-Gutierrez, D. (accepted). The demographic drivers of grief and memory after genocide in Guatemala. Demography.

Reproducible materials

Diego Alburez-Gutierrez

• Structure of the data
• Reproducing results from the paper
  • Table 1
    • Define parameters
    • Define functions
    • Run analysis
  • Table 2
    • Define parameters
    • Define functions
    • Run analysis
  • Table 3
    • Define parameters
  • Run analysis
  • Figure 1
    • Define parameters
    • Define functions
    • Run analysis
    • Faceted plot
  • Figure 2
  • Figure 3
    • Define parameters
    • Define functions
    • Run analysis
    • Plot
  • Figure 4
    • Define parameters
    • Define functions:
    • Run analysis
    • Plot
• Session information

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:

  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):

library(EGM)

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.0.4

## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --

## v ggplot2 3.3.3     v purrr   0.3.4
## v tibble  3.1.0     v dplyr   1.0.5
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1

## Warning: package 'ggplot2' was built under R version 4.0.4

## Warning: package 'tibble' was built under R version 4.0.4

## Warning: package 'readr' was built under R version 4.0.4

## Warning: package 'purrr' was built under R version 4.0.4

## Warning: package 'stringr' was built under R version 4.0.4

## Warning: package 'forcats' was built under R version 4.0.4

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(knitr)

## Warning: package 'knitr' was built under R version 4.0.4

library(kableExtra)

## Warning: package 'kableExtra' was built under R version 4.0.4

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(ggrepel)

## Warning: package 'ggrepel' was built under R version 4.0.4

library(patchwork)

## Warning: package 'patchwork' was built under R version 4.0.4

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).

EGM::final %>% select(-DoD_with_alive) %>% head %>% kable

ind_id	id_cuest	DoD	DoB	alive	died_conflict	sex	sex2	origin_literal	origin_departamento	origin_local	parents	father	mother	spouses	children	cousins	grandparents	grandchildren
1	1_1;3_24;33_30	NA	1954	1	FALSE	2	Female	1	baja verapaz	TRUE	8;9	8	9	2	3;120;121	NA	NA	4;60;122;3507;123
2	1_2;3_25;33_36	1984	1939	0	FALSE	1	Male	coban	alta verapaz	FALSE		NA	NA	1	3;120;121	NA	NA	4;60;122;3507;123
3	1_3;3_26;33_37;34_49;25_3	NA	1976	1	FALSE	2	Female	1	baja verapaz	TRUE	2;1	2	1	5;7	4;60;122;123	24;3;50;25;51;3159;75;2199;468;79;85;2200;80;86;81;87;82;76;88;83;89;84;90	8;9	47
4	1_4;33_41;34_55;3_28	NA	1996	1	FALSE	2	Female	2	baja verapaz	TRUE	7;3	7	3	6	47	4;3507;1538;1539;1540;1541	484;1542;2;1	NA
5	1_5;3_27;33_44;34_58	NA	1974	1	FALSE	1	Male	canchun	baja verapaz	FALSE		NA	NA	3	60;122;123	NA	NA	NA
6	1_6;3_29;33_42;34_56;44_31	NA	1993	1	FALSE	1	Male	2	baja verapaz	TRUE	371;178	371	178	4	47	2720;2712;2721;6;2714;419;420;2713;421;2715;422;2716;423;2717;2718;2719	368;361;172;97	NA

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

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

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

 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())

## `summarise()` has grouped output by 'ag_conflict'. You can override using the `.groups` argument.

# 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())

## `summarise()` has grouped output by 'ag_conflict'. You can override using the `.groups` argument.

  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
  )

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

  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) 

## Warning in kable_styling(.): Please specify format in kable. kableExtra can
## customize either HTML or LaTeX outputs. See https://haozhu233.github.io/
## kableExtra/ for details.

## Warning in add_header_above(., c("", "", `(A)` = 3, `(B)` = 3, `(C) = B/A` =
## 3)): Please specify format in kable. kableExtra can customize either HTML or
## LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.

## Warning in add_header_above(., c(Cohort = 1, `Age in 1982` = 1, `Alive in 1982`
## = 3, : Please specify format in kable. kableExtra can customize either HTML or
## LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.

		All	Female	Male	All	Female	Male	All	Female	Male
1968-1982	0 to 14	544	273	271	163	82	81	.3	.3	.3
1937-1967	15 to 44	332	170	162	147	75	72	.44	.44	.44
Pre-1937	45+	94	43	51	56	22	34	.59	.51	.67
All cohorts	All ages	970	486	484	366	179	187	.38	.37	.39

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.

write.csv(em_df, "Tab1.csv", row.names = F)

Table 2

Define parameters

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

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

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)

## `summarise()` has grouped output by 'sex', 'ag'. You can override using the `.groups` argument.

# 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)

## `summarise()` has grouped output by 'sex'. You can override using the `.groups` argument.

# 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")
  ) 

	F 0-14	F 15-44	F 45+	F Tot	M 0-14	M 15-44	M 45+	M Tot	All_Total
Children	0 (0.1)	0.5 (1.2)	1.8 (1.6)	0.3 (0.9)	0 (0.1)	0.7 (1.3)	1.8 (2.1)	0.3 (1)	0.3 (0.9)
Parents	0.8 (0.7)	0.8 (0.8)	0.6 (0.8)	0.8 (0.7)	0.8 (0.7)	0.9 (0.7)	0.4 (0.7)	0.8 (0.7)	0.8 (0.7)
Siblings	0.6 (1.3)	1.7 (2.5)	1 (1.3)	1 (1.8)	0.8 (1.3)	1.4 (1.7)	0.9 (2.1)	1 (1.5)	1 (1.6)
Spouses	0 (0.1)	0.2 (0.4)	0.5 (0.5)	0.1 (0.3)	0 (0)	0.2 (0.4)	0.3 (0.5)	0.1 (0.3)	0.1 (0.3)
Grandchildren	0 (0)	0.3 (0.8)	2.6 (3.8)	0.3 (1.2)	0 (0)	0.2 (0.8)	2.7 (3.7)	0.2 (1.1)	0.2 (1.2)
Grandparents	0.7 (1)	0.3 (0.7)	0 (0)	0.5 (0.9)	0.8 (1)	0.2 (0.5)	0 (0)	0.5 (0.9)	0.5 (0.9)
Cousins	2 (3.9)	1.4 (3.5)	0 (0)	1.7 (3.7)	2.3 (4)	1.3 (3.6)	0 (0)	1.9 (3.8)	1.8 (3.8)
Nuclear (any)	1.2 (1.7)	2.9 (3.2)	3.5 (3.1)	1.8 (2.5)	1.4 (1.7)	2.9 (2.8)	3.2 (3.3)	1.9 (2.3)	1.9 (2.4)
Extended (Total)	3.8 (5.4)	4.9 (5.3)	6.1 (5.4)	4.3 (5.4)	4.5 (5.5)	4.6 (5.1)	5.8 (6.3)	4.6 (5.4)	4.4 (5.4)

Number of relatives lost by sex and age of genocide survivors at the time of the killings (mean and standard deviation).

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

Table 3

Define parameters

deaths_keep <- 1:5

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

Run analysis

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"
    ) 

(n-1) to n+	0 to 1+	1 to 2+	2 to 3+	3 to 4+	4 to 5+
Nuclear (Any)	.59 [354]	.68 [239]	.79 [188]	.63 [119]	.72 [86]
Children	.11 [68]	.6 [41]	.49 [20]	.35 [7]	.71 [5]
Parents	.42 [250]	.28 [71]
Siblings	.38 [228]	.71 [162]	.67 [108]	.47 [51]	.49 [25]
Spouses	.07 [44]
Extended (Any)	.66 [397]	.86 [342]	.91 [312]	.83 [259]	.84 [217]
Grandchildren	.08 [45]	.69 [31]	.65 [20]	.65 [13]	.69 [9]
Grandparents	.28 [170]	.62 [105]	.3 [32]	.03 [1]
Cousins	.29 [173]	.91 [158]	.78 [124]	.76 [94]	.98 [92]

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.

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

Figure 1

Define parameters

# 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

# 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

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

# 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?

# 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")
  )

## Warning: Removed 2 rows containing missing values (geom_point).

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

## Warning: Removed 2 rows containing missing values (geom_point).

Figure 3

Define parameters

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

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

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

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

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:

# 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

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

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

## [1] "Report created: 2022-05-02 10:48:44"

## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.1252 
## [2] LC_CTYPE=English_United Kingdom.1252   
## [3] LC_MONETARY=English_United Kingdom.1252
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] rionegro_0.1.0   patchwork_1.1.1  ggrepel_0.9.1    kableExtra_1.3.4
##  [5] knitr_1.31       forcats_0.5.1    stringr_1.4.0    dplyr_1.0.5     
##  [9] purrr_0.3.4      readr_1.4.0      tidyr_1.1.3      tibble_3.1.0    
## [13] ggplot2_3.3.3    tidyverse_1.3.0  EGM_0.1.0       
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.7        svglite_2.0.0     lubridate_1.7.10  assertthat_0.2.1 
##  [5] digest_0.6.28     utf8_1.2.1        plyr_1.8.6        R6_2.5.0         
##  [9] cellranger_1.1.0  backports_1.2.1   reprex_1.0.0      evaluate_0.14    
## [13] httr_1.4.2        highr_0.8         pillar_1.5.1      rlang_0.4.10     
## [17] readxl_1.3.1      rstudioapi_0.13   rmarkdown_2.7     labeling_0.4.2   
## [21] webshot_0.5.2     munsell_0.5.0     broom_0.7.5       compiler_4.0.2   
## [25] modelr_0.1.8      xfun_0.21         pkgconfig_2.0.3   systemfonts_1.0.1
## [29] htmltools_0.5.1.1 tidyselect_1.1.0  fansi_0.4.2       viridisLite_0.3.0
## [33] crayon_1.4.1      dbplyr_2.1.0      withr_2.4.1       grid_4.0.2       
## [37] jsonlite_1.7.2    gtable_0.3.0      lifecycle_1.0.0   DBI_1.1.1        
## [41] magrittr_2.0.1    scales_1.1.1      cli_2.3.1         stringi_1.5.3    
## [45] farver_2.1.0      fs_1.5.0          xml2_1.3.2        ellipsis_0.3.1   
## [49] generics_0.1.0    vctrs_0.3.7       tools_4.0.2       glue_1.4.2       
## [53] hms_1.0.0         yaml_2.2.1        colorspace_2.0-0  rvest_1.0.0      
## [57] haven_2.3.1