Nonlinear temporal (conception date) 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(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(wk_ini)", 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(wk_ini)", 1, intercept = FALSE, FUN = c95_exp
    ) +
    labs(y = "LBW odds ratio", x = "Conception date", subtitle = "(b)",
         fill = legend_title, col = legend_title, linetype = legend_title) +
    theme(legend.position = "bottom")

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

gg <- gg_1_1 + gg_1_2 + gg_1_3 +
    plot_layout(guides = "collect", nrow = 1) &
    theme(legend.position = "bottom", plot.margin = margin(3, 7, 0, 0))

# add date labels
rain_date0 <- lubridate::ymd("2002-01-01") # to 2013-12-31
date_at <- 1 + floor(difftime(rain_date0 + lubridate::years(seq(1, 20, 2)),
                              rain_date0, unit = "weeks"))
date_at <- as.numeric(date_at)
date_labels <- 2002 + seq(1, 20, 2)
gg <- gg &
    scale_x_continuous(breaks = as.numeric(date_at), labels = date_labels)

base_size <- 7
gg <- gg &
    theme(
          legend.title = element_text(size = base_size, face = "bold"),
          legend.text = element_text(size = base_size),
          axis.title = element_text(size = base_size),
          axis.text = element_text(size = base_size),
          strip.text = element_text(size = base_size),
    )
ggsave(file.path(path_summarised, "paper-mod-eff-nonlinear-temporal.pdf"), gg ,
       width = 18, height = 6.3, units = "cm")
ggsave(file.path(path_summarised, "paper-mod-eff-nonlinear-temporal.jpg"), gg ,
       width = 18, height = 6.3, units = "cm", dpi = 350)

Visualize nonlinear temporal effects: conception date

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 
#>  50.690   2.947  53.621