Rainfall 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

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"))

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')
baseline <- TRUE

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, baseline = TRUE)

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, baseline = TRUE)

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, baseline = TRUE)

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, baseline = TRUE)

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, baseline = TRUE)

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, model.frame = bwdata_model, baseline = TRUE,
    labs = non_extreme_labs, labs_weeks = non_extreme_labs_weeks)

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, model.frame = bwdata_model, baseline = TRUE,
    labs = non_extreme_labs, labs_weeks = non_extreme_labs_weeks)

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,
    model.frame = bwdata_model, reference = 1, baseline = TRUE,
    labs = non_extreme_labs, labs_weeks = non_extreme_labs_weeks, labs_ci = labs_lbw)

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,
    model.frame = bwdata_model, reference = 1, baseline = TRUE,
    labs = non_extreme_labs, labs_weeks = non_extreme_labs_weeks, labs_ci = labs_lbw)

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,
    model.frame = bwdata_model, reference = 1, baseline = TRUE,
    labs = non_extreme_labs, labs_weeks = non_extreme_labs_weeks, labs_ci = labs_pre)

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 = fun_mu, model.frame = bwdata_model, baseline = TRUE,
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks)

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 = fun_mu, model.frame = bwdata_model, baseline = TRUE,
    labs = sea_dev_labs, labs_weeks = sea_dev_labs_weeks)

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 = 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_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 = 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_sea_dev_pre <- gg_rain(
    model_pre, 'pi', sea_dev_term, grid = grid, data_points = data_points_all[[1]],
    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)

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"
# )
model_labs <- paste0("(", letters[1:15], ")")
labs0 <- labs(x = NULL, y = NULL)
gg <-
    gg_ext_bw_full$significant + labs0 + labs(title = model_labs[1]) +
    gg_non_ext_bw_full$significant + labs0 + labs(title = model_labs[2]) +
    gg_sea_dev_bw_full$significant + labs0 + labs(title = model_labs[3]) +
    gg_ext_bw_growth$significant + labs0 + labs(title = model_labs[4]) +
    gg_non_ext_bw_growth$significant + labs0 + labs(title = model_labs[5]) +
    gg_sea_dev_bw_growth$significant + labs0 + labs(title = model_labs[6]) +
    gg_ext_lbw_full$significant + labs0 +
    labs(y = "Exposure to extreme intense rainfall", title = model_labs[7]) +
    theme(axis.title.y = element_text(hjust = -6)) +
    gg_non_ext_lbw_full$significant + labs0 +
    labs(y = non_extreme_labs[2], title = model_labs[8]) +
    theme(axis.title.y = element_text(hjust = -3)) +
    gg_sea_dev_lbw_full$significant + labs0 +
    labs(y = sea_dev_labs[2], title = model_labs[9]) +
    theme(axis.title.y = element_text(hjust = -8.5)) +
    gg_ext_lbw_growth$significant + labs0 + labs(title = model_labs[10]) +
    gg_non_ext_lbw_growth$significant + labs0 + labs(title = model_labs[11]) +
    gg_sea_dev_lbw_growth$significant + labs0 + labs(title = model_labs[12]) +
    gg_ext_pre$significant +
    labs(x = day2day::wrapit("Exposure to extreme deficient rainfall", 30),
         y = "Exposure to extreme intense rainfall", title = model_labs[13]) +
    theme(axis.title.y = element_text(hjust = -6)) +
    gg_non_ext_pre$significant +
    labs(x = day2day::wrapit(non_extreme_labs[1], 30),
         y = non_extreme_labs[2], title = model_labs[14]) +
    theme(axis.title.y = element_text(hjust = -3)) +
    gg_sea_dev_pre$significant +
    labs(x = day2day::wrapit(sea_dev_labs[1], 30),
         y = sea_dev_labs[2], title = model_labs[15]) +
    theme(axis.title.y = element_text(hjust = -8.5)) +
    plot_layout(nrow = 5, guides = "collect")

base_size <- 7
gg <- gg &
    theme(
          strip.background = element_blank(),
          strip.text.x = element_blank(),
          plot.title = element_text(size = base_size),
          axis.text = element_text(size = base_size),
          axis.title = element_text(size = base_size),
          plot.margin = unit(c(0, 0.25, 0.05, 0.3), "cm")
          # axis.title.y = element_text(hjust = 0)
    )
ggsave(file.path(path_summarised, "paper-mod-eff-rain-baseline-significant.pdf"), gg,
       width = 18, height = 21, units = "cm")

ggsave(file.path(path_summarised, "paper-mod-eff-rain-baseline-significant.jpg"), gg,
       width = 18, height = 21, units = "cm", dpi = 350)
# 20 27

Visualize significant effects for BW, LBW and PTB models

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

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

model_labs <- paste0("(", letters[1:5], ")")
rain_effect_labs <- c("Extreme 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))

base_size <- 7

ggp_ci <- function(data, trans = "identity", base_size = 7) {
    gg_ci <- ggplot(data, aes(description, mean)) +
        geom_errorbar(aes(ymin = lower, ymax = upper, color = rain_effect),
                      position = position_dodge(0.3), size = rel(0.4), width = 0.25) +
        geom_point(aes(color = rain_effect),
                   position = position_dodge(0.3), size = rel(1)) +
        # 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(base_size = base_size, base_family = "Helvetica") +
        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",
          legend.title = element_text(size = base_size, face = "bold"),
          legend.text = element_text(size = base_size)
          ) +
    plot_layout(ncol = 2) &
    theme(
          strip.background = element_blank(),
          strip.text.x = element_text(hjust = -0.01, size = base_size),
          axis.text = element_text(size = base_size)
    )

ggsave(file.path(path_summarised, "paper-mod-eff-rain-baseline-ci.pdf"), gg_ci,
       width = 18, height = 18, units = "cm")

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

print(gg_ci, 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 
#> 168.035   4.179 172.275