Nonlinear effects of maternal age and hydrological seasonality 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(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)
library(patchwork)

#> 
#> Attaching package: 'patchwork'

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

library(ggplot2)

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", "56-bamlss-vis-gg.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)

Create figure

types <- factor(c("No", "Yes"))
legend_title <- "Controlling for gestational age?"
guides_none <- guides(fill = "none", col = "none", linetype = "none")

gg_1_1 <- gg_smooths(
    list(model_bw_full, model_bw_growth), bwdata_model,
    "mu", "s(age)", 0, intercept = FALSE
    ) +
    labs(y = "Mean change (grams)", subtitle = "(a)") +
    guides_none

#> Note: Using an external vector in selections is ambiguous.
#> ℹ Use `all_of(var)` instead of `var` to silence this message.
#> ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
#> This message is displayed once per session.

gg_1_2 <- gg_smooths(
    list(model_lbw_full, model_lbw_growth), bwdata_model,
    "pi", "s(age)", 1, intercept = FALSE, FUN = c95_exp
    ) +
    labs(y = "LBW odds ratio", x = "Age", subtitle = "(b)") +
    guides_none

gg_1_3 <- gg_smooths(
    list(model_pre), bwdata_model,
    "pi", "s(age)", 1, types = types[1], intercept = FALSE, FUN = c95_exp
    ) +
    labs(y = "PTB odds ratio", subtitle = "(c)") +
    guides_none

gg_2_1 <- gg_smooths_river(
    list(model_bw_full, model_bw_growth), bwdata_model,
    "mu", "s(rivwk_conception)", 0, intercept = FALSE
    ) +
    labs(y = "Mean change (grams)", subtitle = "(d)") +
    guides_none

gg_2_2 <- gg_smooths_river(
    list(model_lbw_full, model_lbw_growth), bwdata_model,
    "pi", "s(rivwk_conception)", 1, intercept = FALSE, FUN = c95_exp
    ) +
    labs(y = "LBW odds ratio", x = "Seasonal river level index", subtitle = "(e)",
         fill = legend_title, col = legend_title, linetype = legend_title) +
    theme(legend.position = "bottom")

gg_2_3 <- gg_smooths_river(
    list(model_pre), bwdata_model,
    "pi", "s(rivwk_conception)", 1, types = types[1], intercept = FALSE, FUN = c95_exp
    ) +
    labs(y = "PTB odds ratio", subtitle = "(f)") +
    guides_none

base_size <- 7
gg <- gg_1_1 + gg_1_2 + gg_1_3 + gg_2_1  + gg_2_2 + gg_2_3 +
    plot_layout(guides = "collect", nrow = 2) &
    theme(legend.position = "bottom", plot.margin = margin(3, 5, 0, 0),
          legend.title = element_text(size = base_size, face = "bold"),
          legend.text = element_text(size = base_size),
          axis.text = element_text(size = base_size),
          axis.title = element_text(size = base_size),
          plot.subtitle = element_text(size = base_size),
          axis.text.x.top = element_text(colour = "steelblue2"))

ggsave(file.path(path_summarised, "paper-mod-eff-nonlinear-age-and-seasonality.pdf"), gg,
       width = 18, height = 12, units = "cm")

Visualize nonlinear effects of maternal age and hydrological seasonality

print(gg, 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 
#>  53.082   2.897  55.956