# Rainfall effects on BW, LBW and PTB without baselines

In this script, we create figures to summarise the rainfall variability effects of the BW, LBW, and PTB models without using pre-defined baselines. This visualization might not be adequate given that the implicit baselines do not represent mother without or low exposure. The results can be seen at:

- Section Visualize significant effects for BW, LBW and PTB
models, which is an alternative visualization of the
*Extended Data Figure 7*of our paper. - Section Visualize reference points effects (with credible intervals) for BW, LBW and
PTB models, which is an alternative visualization of the
*Figure 4*of our paper.

## 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-nobaseline-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-nobaseline-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
#> 143.415 3.882 147.355
```