Counts of newborn’s prenatal exposure

In this script, we compute the counts of newborn’s prenatal exposure to (a) seasonal deviations in rainfall, (b) non-extreme rainfall events, and (c) extreme rainfall events. The resulting figure can be seen at section Visualize counts of newborn’s prenatal exposure. This output corresponds to Extended Data Fig. 3 of our paper.

This analysis helps to define common groups of exposure to be used as baselines.

Load packages, read data and source custom scripts

rm(list = ls())
library(ggplot2)
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(patchwork)

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

source(file.path(path_proj, "src", "60-count2d.R"))

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

# quantile(abs(bwdata_model$neg_ext_mbsi_mean_8wk) + abs(bwdata_model$pos_ext_mbsi_mean_8wk),
#          probs = seq(0, 1, 0.1))

Code to compute counts by levels of exposure

text_color <- "red"
point_color <- "red"
bins <- 50

# extremes
situations_ext <- c("None", "Intense and deficient", "Intense", "Moderate intense",
                    "Deficient")
data_points_ext <- data.frame(
    x = c(-0.0, -0.51, -0.02, -0.04, -1.03),
    y = c(0.0, 0.37, 0.78, 0.5, 0.03),
    label = c("A", "B", "C", "E", "D"),
    description = factor(situations_ext, levels = situations_ext),
    vjust = c(0, 0.3, 0.3, 0.3, -1),
    hjust = c(1.5, 1.5, 1.5, 1.5, 0.5)
)
count_ext <- exposure2d_count(
    bwdata_model,
    c("neg_ext_mbsi_mean_8wk", "pos_ext_mbsi_mean_8wk"),
    data_points = data_points_ext, bins,
    point_color = point_color, text_color = text_color
)
count_ext <- count_ext +
    labs(x = 'Extreme deficient rainfall',
         y = 'Extreme intense rainfall',
         fill = 'Counts',
         subtitle = '(c)')

# non-extremes
situations <- c("None", "Intense and deficient", "Intense", "Deficient")
data_points_non_ext <- data.frame(
    x = c(-0.25, -0.65, -0.02, -1.6),
    y = c(0.4, 0.65, 1.2, 0.03),
    label = c("A", "B", "C", "D"),
    description = factor(situations, levels = situations),
    vjust = c(0.3, 0.3, 0.3, -1),
    hjust = c(1.5, 1.5, 1.5, 0.5)
)
count_non_ext <- exposure2d_count(
    bwdata_model,
    c("neg_mckee_mean_8wk", "pos_mckee_mean_8wk"),
    data_points = data_points_non_ext, bins,
    point_color = point_color, text_color = text_color
)
count_non_ext <- count_non_ext +
    labs(x = 'Non-extreme deficient rainfall',
         y = 'Non-extreme intense rainfall',
         fill = 'Counts',
         subtitle = '(b)')

# seasonal deviation
data_points_sea_dev <- data.frame(
    x = c(-0.37, -0.55, -0.23, -0.8),
    y = c(0.4, 0.55, 0.7, 0.12),
    label = c("A", "B", "C", "D"),
    description = factor(situations, levels = situations),
    vjust = c(0.3, 0.3, 0.3, 0.3),
    hjust = c(1.5, 1.5, 1.5, 1.5)
)
count_deviation <- exposure2d_count(
    bwdata_model,
    c("neg_mbsi_mean_1wk", "pos_mbsi_mean_1wk"),
    data_points = data_points_sea_dev, bins,
    point_color = point_color, text_color = text_color
)
count_deviation <- count_deviation +
    labs(x = 'Negative deviation from seasonality',
         y = 'Positive deviation from seasonality',
         fill = 'Counts',
         subtitle = '(a)')

Code to create figure of counts of newborn’s prenatal exposure

base_size <- 7
theme <- theme_bw(base_size = base_size, base_family = "Helvetica") +
    theme(legend.key.width = unit(1, "lines"),
          legend.key.height = unit(0.5, "lines"),
          legend.title = element_text(size = base_size, vjust = 1, face = "bold"),
          legend.text = element_text(size = base_size),
          legend.position = "bottom",
          axis.text = element_text(size = base_size),
          axis.title = element_text(size = base_size),
          strip.text = element_text(size = base_size),
          plot.margin = unit(c(0.1, 0.3, 0, 0.2), "lines"),
    )

gg <- (count_deviation | count_non_ext | count_ext) &
    theme

ggsave(file.path(path_summarised, "paper-data-exposure-counts.pdf"), gg,
       width = 18, height = 6.2, units = "cm")
ggsave(file.path(path_summarised, "paper-data-exposure-counts.jpg"), gg,
       width = 18, height = 6.2, units = "cm", dpi = 350)

Visualize counts of newborn’s prenatal exposure

gg + plot_layout(ncol = 1, nrow = 3)

Save reference points

saveRDS(list(data_points_sea_dev, data_points_non_ext, data_points_ext),
        file.path(path_processed, "exposure-reference-points.rds"))

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  23.763   0.813  24.559