Visualization of the effects of extreme events using empirical data

Load packages, read data and source custom scripts

Paths are defined relative to the git repository location.

rm(list = ls())
library(dplyr)

#> 
#> Attaching package: 'dplyr'

#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag

#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

library(ggplot2)

path_proj <- day2day::git_path()
path_data <- file.path(path_proj, "data")
path_processed <- file.path(path_data, "processed")
path_out_sum_data <- file.path(path_proj, "data", "summarised")

bwdata_model <- fst::read_fst(file.path(path_processed, "bwdata_41_model.fst"))

Organize data to make visualization

bwdata_plot <- bwdata_model %>%
    mutate(group1 = cut(neg_ext_mbsi_mean_8wk, c(0.1, -0.154, -0.309, -0.463,
                                                 -0.926, -2)),
           group2 = cut(pos_ext_mbsi_mean_8wk, c(-0.1, 0.119, 0.357, 0.714, 1)),
           group3 = consult_num
    ) %>%
    within({
        levels(group1) <- paste0(c("(High)", "(Moderate)", "(Intermediate)", "(Low)", "(None)"),
        "\n exposure to extremely \n deficient events")
        levels(group2) <- paste0(rev(c("(High)", "(Moderate)", "(Intermediate)", "(None)")),
        "\n exposure to extremely \n intense events")
    })

bwdata_plot <- bwdata_plot %>%
    group_by(group1, group2, group3) %>%
    summarise(n = n(), bw_mean = mean(born_weight), bw_sd = sd(born_weight),
              lbw_mean = mean(lbw), lbw_sd = sqrt(lbw_mean * (1 - lbw_mean))) %>%
    mutate(bw_se = bw_sd / sqrt(n),
           bw_lower = bw_mean - 1.96 * bw_se,
           bw_upper = bw_mean + 1.96 * bw_se,
           lbw_se = lbw_sd / sqrt(n),
           lbw_lower = lbw_mean - 1.96 * lbw_se,
           lbw_upper = lbw_mean + 1.96 * lbw_se
           ) %>%
    tidyr::gather(varname, value, bw_mean:lbw_upper) %>%
    tidyr::extract(varname, c("type", "stat"), "([[:alnum:]]+)_([[:alnum:]]+)") %>%
    tidyr::spread(stat, value) %>%
    mutate(type = factor(type, c("lbw", "bw"),
                         c("Proportion of low birthweight", "Mean birthweight")))

#> `summarise()` regrouping output by 'group1', 'group2' (override with `.groups` argument)

bwdata_bw <- bwdata_plot %>%
    filter(type == "Mean birthweight")

Create and save figure

base_size <- 7
gg1 <- ggplot(subset(bwdata_bw, !is.na(lower))) +
    geom_errorbar(aes(x = group3, ymin = lower, ymax = upper, color = group3),
                  show.legend = FALSE, width = 0.3) +
    geom_point(aes(x = group3, y = mean), col = 1, size = rel(0.4)) +
    facet_grid(group1 ~ group2) +
    geom_hline(yintercept = 3220, col = 1, linetype = 2, size = rel(0.3)) +
    labs(color = NULL,
         x = "Number of antenatal consultations",
         y = "Birth-weight (grams)") +
    theme_bw(base_size = base_size, base_family = "Helvetica") +
    theme(
          panel.grid.minor = element_blank(),
          axis.title = element_text(size = base_size),
          axis.text = element_text(size = base_size),
          strip.text = element_text(size = base_size),
          strip.background = element_rect(fill = "gray96")
    )

ggsave(file.path(path_out_sum_data, "paper-data-bw-per-exposure.pdf"),
       gg1, width = 18, height = 18, units = "cm")
ggsave(file.path(path_out_sum_data, "paper-data-bw-per-exposure.jpg"),
       gg1, width = 18, height = 18, units = "cm", dpi = 350)

Visualize birth weight (grams) and under different levels of exposure

print(gg1, vp = grid::viewport(gp = grid::gpar(cex = 1.15)))

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  17.547   0.563  18.090