Extreme events effects on BW, LBW and PTB


In this script, we create figures to summarise the extreme events effects of the BW, LBW, and PTB models. These figures are related to Figure 4 and Extended Data Figure 7 of our paper.

Load packages, read data and source custom scripts

rm(list = ls())
library(bamlss)
#> Loading required package: coda
#> Loading required package: colorspace
#> Loading required package: mgcv
#> Loading required package: nlme
#> This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
#> 
#> Attaching package: 'bamlss'
#> The following object is masked from 'package:mgcv':
#> 
#>     smooth.construct
library(gamlss.dist)
#> Loading required package: MASS
library(ggplot2)
library(patchwork)
#> 
#> Attaching package: 'patchwork'
#> The following object is masked from 'package:MASS':
#> 
#>     area
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", "51-bamlss.R"))
source(file.path(path_proj, "src", "57-bamlss-vis-paper.R"))
source(file.path(path_proj, "src", "53-vis-bw-lbw.R"))

bwdata_file <- file.path(path_processed, "bwdata_41_model.fst")
path_model_bw_full <- file.path(path_modelled, "bw-10-full-re-t-burned.rds")
path_model_bw_growth <- file.path(path_modelled, "bw-20-growth-re-t-burned.rds")
path_model_lbw_full <- file.path(path_modelled, "lbw-10-full-re-burned.rds")
path_model_lbw_growth <- file.path(path_modelled, "lbw-20-growth-re-burned.rds")
path_model_pre <- file.path(path_modelled, "pre-11-full-re-burned.rds")

bwdata_model <- fst::read_fst(bwdata_file)
model_bw_full <- readRDS(path_model_bw_full)
model_bw_growth <- readRDS(path_model_bw_growth)
model_lbw_full <- readRDS(path_model_lbw_full)
model_lbw_growth <- readRDS(path_model_lbw_growth)
model_pre <- readRDS(path_model_pre)

data_points <- readRDS(file.path(path_processed, "exposure-reference-points.rds"))[[3]]

Create set of figures for BW, LBW and PTB

fun_mu <- function (x) c95_expression(x, prefix = paste0('(', letters[1:3], ') '))
fun_pi <- function (x) c95_expression(exp(x), 'odds')
labs_lbw <- c("Prenatal exposure to extreme rainfall events", "Odds ratio of low birthweight")
labs_pre <- c("Prenatal exposure to extreme rainfall events", "Odds ratio of prematurity")
baseline <- TRUE
grid <- 50
term <- 's(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk)'

gg_bw_full <- gg_rain(
    model_bw_full, 'mu', term, grid = grid, data_points = data_points,
    binwidth = 20, FUN = fun_mu, model.frame = bwdata_model, baseline = TRUE)

gg_bw_growth <- gg_rain(
    model_bw_growth, 'mu', term, grid = grid, data_points = data_points,
    binwidth = 5, FUN = fun_mu, model.frame = bwdata_model, baseline = TRUE)

data_points_bw <- data_points
data_points_bw[data_points_bw$label == "E", "y"] <- 0.4

gg_lbw_full <- gg_rain(
    model_lbw_full, 'pi', term, grid = grid, data_points = data_points_bw,
    binwidth = 0.5, FUN = fun_pi, trans = 'log', rev = TRUE,
    model.frame = bwdata_model, reference = 1, baseline = TRUE,
    labs_ci = labs_lbw)

gg_lbw_growth <- gg_rain(
    model_lbw_growth, 'pi', term, grid = grid, data_points = data_points_bw,
    binwidth = 0.02, FUN = fun_pi, trans = 'log', rev = TRUE,
    labs_ci = labs_lbw,
    model.frame = bwdata_model, reference = 1, baseline = TRUE)

gg_pre <- gg_rain(
    model_pre, 'pi', term, grid = grid, data_points = data_points,
    binwidth = 2, FUN = fun_pi, trans = 'log', rev = TRUE,
    labs_ci = labs_pre,
    model.frame = bwdata_model, reference = 1, baseline = TRUE)

Main effects

gg <-
    (gg_bw_full$quantiles +
     labs(x = NULL, y = NULL, subtitle = "1) Full birthweight model")) +
    (gg_bw_growth$quantiles +
     labs(x = NULL, y = NULL, subtitle = "2) Growth birthweight model")) +
    (gg_lbw_full$quantiles +
     labs(x = NULL, subtitle = "3) Full low birthweight model") +
     theme(axis.title.y = element_text(hjust = 0))) +
    (gg_lbw_growth$quantiles +
     labs(x = NULL, y = NULL, subtitle = "4) Growth low birthweight model")) +
    (gg_pre$quantiles +
     labs(y = NULL, subtitle = "5) Prematurity model")) +
    plot_layout(ncol = 1)
ggsave(file.path(path_summarised, "mod-eff-rain-extremes.pdf"), gg ,
       width = 20, height = 32, units = "cm")

print(gg)

Main significant effects

gg <-
    (gg_bw_full$significant +
     labs(x = NULL, subtitle = "1) Full birthweight model")) +
    (gg_bw_growth$significant +
     labs(x = NULL, y = NULL, subtitle = "2) Growth birthweight model")) +
    (gg_lbw_full$significant +
     labs(x = NULL, subtitle = "3) Full low birthweight model")) +
    (gg_lbw_growth$significant +
     labs(x = NULL, y = NULL, subtitle = "4) Growth low birthweight model")) +
    (gg_pre$significant + labs(subtitle = "5) Prematurity model")) +
    plot_layout(ncol = 2, nrow = 3)
ggsave(file.path(path_summarised, "mod-eff-rain-extremes-significant.pdf"), gg ,
       width = 21, height = 25, units = "cm")

print(gg)

Credible intervals for reference points

gg <-
    (gg_bw_full$ci +
     labs(x = NULL, subtitle = "1) Full birthweight model")) +
    (gg_bw_growth$ci +
     labs(x = NULL, y = NULL, subtitle = "2) Growth birthweight model")) +
    (gg_lbw_full$ci +
     labs(x = NULL, subtitle = "3) Full low birthweight model")) +
    (gg_lbw_growth$ci +
     labs(x = NULL, y = NULL, subtitle = "4) Growth low birthweight model")) +
    (gg_pre$ci + labs(subtitle = "5) Prematurity model")) +
    plot_layout(ncol = 2, nrow = 3)
ggsave(file.path(path_summarised, "mod-eff-rain-extremes-ci.pdf"),
       gg, width = 27, height = 21, units = "cm")

print(gg)

Main effects (weeks approximation)

gg <-
    (gg_bw_full$quantiles_weeks +
     labs(x = NULL, y = NULL, subtitle = "1) Full birthweight model")) +
    (gg_bw_growth$quantiles_weeks +
     labs(x = NULL, y = NULL, subtitle = "2) Growth birthweight model")) +
    (gg_lbw_full$quantiles_weeks +
     labs(x = NULL, subtitle = "3) Full low birthweight model") +
     theme(axis.title.y = element_text(hjust = 0))) +
    (gg_lbw_growth$quantiles_weeks +
     labs(x = NULL, y = NULL, subtitle = "4) Growth low birthweight model")) +
    (gg_pre$quantiles_weeks +
     labs(y = NULL, subtitle = "5) Prematurity model")) +
    plot_layout(ncol = 1)
ggsave(file.path(path_summarised, "mod-eff-rain-extremes-weeks.pdf"), gg ,
       width = 20, height = 32, units = "cm")

print(gg)

Main significant effects (weeks approximation)

gg <-
    (gg_bw_full$significant_weeks +
     labs(x = NULL, y = NULL, subtitle = "1) Full birthweight model")) +
    (gg_bw_growth$significant_weeks +
     labs(x = NULL, y = NULL, subtitle = "2) Growth birthweight model")) +
    (gg_lbw_full$significant_weeks +
     labs(x = NULL, subtitle = "3) Full low birthweight model") +
     theme(axis.title.y = element_text(hjust = 0))) +
    (gg_lbw_growth$significant_weeks +
     labs(x = NULL, y = NULL, subtitle = "4) Growth low birthweight model")) +
    (gg_pre$significant_weeks +
     labs(y = NULL, subtitle = "5) Prematurity model")) +
    plot_layout(ncol = 2, nrow = 3)
ggsave(file.path(path_summarised, "mod-eff-rain-extremes-significant-weeks.pdf"),
       gg, width = 21, height = 25, units = "cm")

print(gg)

Time to execute the task

Only useful when executed with Rscript.

proc.time()
#>    user  system elapsed 
#> 150.405   3.798 154.107