Seasonal deviations effects on BW, LBW and PTB

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"))[[1]]

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_mbsi_mean_1wk,pos_mbsi_mean_1wk)"
sea_dev_labs <- c('Negative deviation from seasonality', 'Positive deviation from seasonality')
sea_dev_labs_weeks <- c('Weeks of negative deviation from seasonality',
                        'Weeks of positive deviation from seasonality')

gg_bw_full <- gg_rain(
    model_bw_full, 'mu', term, grid = grid, data_points = data_points,
    binwidth = 5, skip = 0, FUN = fun_mu, model.frame = bwdata_model, baseline = baseline,
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks)

gg_bw_growth <- gg_rain(
    model_bw_growth, 'mu', term, grid = grid, data_points = data_points,
    binwidth = 5, skip = 0, FUN = fun_mu, model.frame = bwdata_model, baseline = baseline,
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks)

gg_lbw_full <- gg_rain(
    model_lbw_full, 'pi', term, grid = grid, data_points = data_points,
    skip = 0, binwidth = 0.05, FUN = fun_pi, trans = 'log', rev = TRUE,
    model.frame = bwdata_model, reference = 1, baseline = TRUE,
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks, labs_ci = labs_lbw)

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

gg_pre <- gg_rain(
    model_pre, 'pi', term, grid = grid, data_points = data_points,
    skip = 0, binwidth = 0.5, FUN = fun_pi, trans = 'log', rev = TRUE,
    model.frame = bwdata_model, reference = 1, baseline = TRUE,
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks, labs_ci = labs_pre)

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-seasonal-deviations.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-seasonal-deviations-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-seasonal-deviations-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-seasonal-deviations-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-seasonal-deviations-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 
#> 155.393   3.845 159.513