Rainfall effects on BW, LBW and PTB without baselines

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

library(dplyr)

#> 
#> Attaching package: 'dplyr'

#> The following object is masked from 'package:MASS':
#> 
#>     select

#> The following object is masked from 'package:bamlss':
#> 
#>     n

#> The following object is masked from 'package:nlme':
#> 
#>     collapse

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

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

library(purrr)

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_all <- readRDS(file.path(path_processed, "exposure-reference-points.rds"))

data_points_all[[3]][data_points_all[[3]]$label == "A", c("x", "y")] <- c(0.0, 0.0)

Extreme effects

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")
fun_mu <- function (x) c95_expression(x, prefix = paste0('(', letters[1:3], ') '))
fun_pi <- function (x) c95_expression(exp(x), 'odds')

grid <- 50

extreme_term <- 's(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk)'

gg_ext_bw_full <- gg_rain(
    model_bw_full, 'mu', extreme_term, grid = grid, data_points = data_points_all[[3]],
    binwidth = 20, FUN = fun_mu, model.frame = bwdata_model)

gg_ext_bw_growth <- gg_rain(
    model_bw_growth, 'mu', extreme_term, grid = grid, data_points = data_points_all[[3]],
    binwidth = 5, FUN = fun_mu, model.frame = bwdata_model)

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

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

gg_ext_lbw_growth <- gg_rain(
    model_lbw_growth, 'pi', extreme_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)

gg_ext_pre <- gg_rain(
    model_pre, 'pi', extreme_term, grid = grid, data_points = data_points_all[[3]],
    binwidth = 2, FUN = fun_pi, trans = 'log', rev = TRUE,
    labs_ci = labs_pre,
    model.frame = bwdata_model, reference = 1)

Non-extreme effects

non_extreme_term <- "s(neg_mckee_mean_8wk,pos_mckee_mean_8wk)"
non_extreme_labs <- c('Exposure to non-extreme deficient rainfall',
                      'Exposure to non-extreme intense rainfall')
non_extreme_labs_weeks <- c('Weeks exposed to non-extreme deficient rainfall',
                            'Weeks exposed to non-extreme intense rainfall')

gg_non_ext_bw_full <- gg_rain(
    model_bw_full, 'mu', non_extreme_term, grid = grid, data_points = data_points_all[[2]],
    binwidth = 5, FUN = fun_mu,
    labs = non_extreme_labs, labs_weeks = non_extreme_labs_weeks,
    model.frame = bwdata_model)

gg_non_ext_bw_growth <- gg_rain(
    model_bw_growth, 'mu', non_extreme_term, grid = grid, data_points = data_points_all[[2]],
    binwidth = 5, FUN = fun_mu,
    labs = non_extreme_labs, labs_weeks = non_extreme_labs_weeks,
    model.frame = bwdata_model)

gg_non_ext_lbw_full <- gg_rain(
    model_lbw_full, 'pi', non_extreme_term, grid = grid, data_points = data_points_all[[2]],
    binwidth = 0.1, FUN = fun_pi, trans = 'log', rev = TRUE,
    labs = non_extreme_labs, labs_weeks = non_extreme_labs_weeks, labs_ci = labs_lbw,
    model.frame = bwdata_model, reference = 1)

gg_non_ext_lbw_growth <- gg_rain(
    model_lbw_growth, 'pi', non_extreme_term, grid = grid, data_points = data_points_all[[2]],
    binwidth = 0.02, FUN = fun_pi, trans = 'log', rev = TRUE,
    labs = non_extreme_labs, labs_weeks = non_extreme_labs_weeks, labs_ci = labs_lbw,
    model.frame = bwdata_model, reference = 1)

gg_non_ext_pre <- gg_rain(
    model_pre, 'pi', non_extreme_term, grid = grid, data_points = data_points_all[[2]],
    binwidth = 0.5, FUN = fun_pi, trans = 'log', rev = TRUE,
    labs = non_extreme_labs, labs_weeks = non_extreme_labs_weeks, labs_ci = labs_pre,
    model.frame = bwdata_model, reference = 1)

Deviations from seasonality effects

sea_dev_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_sea_dev_bw_full <- gg_rain(
    model_bw_full, 'mu', sea_dev_term, grid = grid, data_points = data_points_all[[1]],
    binwidth = 5,
    FUN = function (x) c95_expression(x, prefix = paste0('(', letters[1:3], ') ')),
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks,
    model.frame = bwdata_model)

gg_sea_dev_bw_growth <- gg_rain(
    model_bw_growth, 'mu', sea_dev_term, grid = grid, data_points = data_points_all[[1]],
    binwidth = 5,
    FUN = function (x) c95_expression(x, prefix = paste0('(', letters[1:3], ') ')),
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks,
    model.frame = bwdata_model)

gg_sea_dev_lbw_full <- gg_rain(
    model_lbw_full, 'pi', sea_dev_term, grid = grid, data_points = data_points_all[[1]],
    binwidth = 0.05,
    FUN = function (x) c95_expression(exp(x), 'odds'), trans = 'log', rev = TRUE,
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks, labs_ci = labs_lbw,
    model.frame = bwdata_model, reference = 1)

gg_sea_dev_lbw_growth <- gg_rain(
    model_lbw_growth, 'pi', sea_dev_term, grid = grid, data_points = data_points_all[[1]],
    binwidth = 0.02,
    FUN = function (x) c95_expression(exp(x), 'odds'), trans = 'log', rev = TRUE,
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks, labs_ci = labs_lbw,
    model.frame = bwdata_model, reference = 1)

gg_sea_dev_pre <- gg_rain(
    model_pre, 'pi', sea_dev_term, grid = grid, data_points = data_points_all[[1]],
    binwidth = 0.5,
    FUN = function (x) c95_expression(exp(x), 'odds'), trans = 'log', rev = TRUE,
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks, labs_ci = labs_pre,
    model.frame = bwdata_model, reference = 1)

Create figure of significant effects for BW, LBW and PTB models

model_labs <- c(
    "1) Full birthweight model",
    "2) Growth birthweight model",
    "3) Full low birthweight model",
    "4) Growth low birthweight model",
    "5) Prematurity model"
)
gg <-
    gg_ext_bw_full$significant + labs(title = model_labs[1]) +
    gg_non_ext_bw_full$significant +
    gg_sea_dev_bw_full$significant +
    gg_ext_bw_growth$significant + labs(title = model_labs[2]) +
    gg_non_ext_bw_growth$significant +
    gg_sea_dev_bw_growth$significant +
    gg_ext_lbw_full$significant + labs(title = model_labs[3]) +
    gg_non_ext_lbw_full$significant +
    gg_sea_dev_lbw_full$significant +
    gg_ext_lbw_growth$significant + labs(title = model_labs[4]) +
    gg_non_ext_lbw_growth$significant +
    gg_sea_dev_lbw_growth$significant +
    gg_ext_pre$significant + labs(title = model_labs[5]) +
    gg_non_ext_pre$significant +
    gg_sea_dev_pre$significant +
    plot_layout(nrow = 5)
gg <- gg &
    theme(
          strip.background = element_blank(),
          strip.text.x = element_blank()
    )
ggsave(file.path(path_summarised, "mod-eff-rain-main-significant.pdf"), gg ,
       width = 25, height = 36, units = "cm")

Visualize significant effects for BW, LBW and PTB models

gg

Create reference points effects (with credible intervals) for BW, LBW and PTB models

rain_effect_labs <- c("Extremes events", "Non-extreme events", "Seasonal deviation")
ci_data <- list(gg_ext_bw_full, gg_non_ext_bw_full, gg_sea_dev_bw_full,
                gg_ext_bw_growth, gg_non_ext_bw_growth, gg_sea_dev_bw_growth,
                gg_ext_lbw_full, gg_non_ext_lbw_full, gg_sea_dev_lbw_full,
                gg_ext_lbw_growth, gg_non_ext_lbw_growth, gg_sea_dev_lbw_growth,
                gg_ext_pre, gg_non_ext_pre, gg_sea_dev_pre) %>%
    map(~ .$ci_data[, - c(1:2)]) %>%
    map2(rep(1:5, each = 3), ~ mutate(.x, model = .y)) %>%
    map2(rep(1:3, 5), ~ mutate(.x, rain_effect = .y)) %>%
    bind_rows() %>%
    mutate(
           model = factor(model, labels = model_labs),
           rain_effect = factor(rain_effect, labels = rain_effect_labs)
    )

intercept_data <- expand.grid(rain_effect = rain_effect_labs, model = model_labs) %>%
    mutate(yintercept = rep(c(0, 0, 1, 1, 1), each = 3))

ggp_ci <- function(data, trans = "identity") {
    gg_ci <- ggplot(data, aes(description, mean)) +
        geom_errorbar(aes(ymin = lower, ymax = upper, color = rain_effect),
                      position = position_dodge(0.3), width = 0.2,
                      size = rel(0.4), width = 0.25) +
        geom_point(aes(color = rain_effect),
                   position = position_dodge(0.3), size = rel(2)) +
        # scale_y_continuous(labels = function(x) round_odds(exp(x))) +
        geom_hline(yintercept = 0, lty = 2, col = 1) +
        labs(x = NULL, y = NULL, color = "Rain effect") +
        facet_wrap(~ model) +
        theme_bw(10) +
        theme(legend.position = "none")
        if (trans == 'log') {
            gg_ci <- gg_ci + scale_y_continuous(labels = function(x) round_odds(exp(x)))
        }
    return(gg_ci)
}

gg_ci <- ggp_ci(subset(ci_data, model == model_labs[1])) +
    ggp_ci(subset(ci_data, model == model_labs[2])) +
    ggp_ci(subset(ci_data, model == model_labs[3]), "log") +
    ggp_ci(subset(ci_data, model == model_labs[4]), "log") +
    ggp_ci(subset(ci_data, model == model_labs[5]), "log") +
    theme(legend.position = "bottom") +
    plot_layout(ncol = 2)

#> Warning: Duplicated aesthetics after name standardisation: width

#> Warning: Duplicated aesthetics after name standardisation: width

#> Warning: Duplicated aesthetics after name standardisation: width

#> Warning: Duplicated aesthetics after name standardisation: width

#> Warning: Duplicated aesthetics after name standardisation: width

ggsave(file.path(path_summarised, "mod-eff-rain-main-ci.pdf"), gg_ci,
       width = 23, height = 20, units = "cm")

Visualize reference points effects (with credible intervals) for BW, LBW and PTB models

gg_ci

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#> 142.793   3.515 146.427