Municipality random effects on BW, LBW and PTB


In this script, we create figures to show the municipality random effects on BW, LBW, and PTB models. The images can be seen at:

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(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
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", "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)

ama <- st_read(file.path(path_processed, "ama_05_62.gpkg")) %>%
    dplyr::select(mun_code, inc_study, mun_name)
#> Reading layer `ama_05_62' from data source `/home/rstudio/Documents/Projects/data-amazonia/70-studies/birthweight/data/processed/ama_05_62.gpkg' using driver `GPKG'
#> Simple feature collection with 62 features and 373 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 1573261 ymin: 8909613 xmax: 3542869 ymax: 10248330
#> projected CRS:  unnamed

Obtain municipality random effects

muni_data <- tibble::tibble(res_muni = unique(bwdata_model$res_muni))

muni_data_bw <- gg_munis(list(model_bw_full, model_bw_growth), muni_data,
                     "mu", "s(res_muni)", FUN = c95, intercept = FALSE, sf = ama)

muni_data_lbw <- gg_munis(list(model_lbw_full, model_lbw_growth), muni_data,
                     "pi", "s(res_muni)", FUN = c95_exp, intercept = FALSE, sf = ama)

muni_data_pre <- gg_munis(list(model_pre), muni_data, types = factor("No"),
                     "pi", "s(res_muni)", FUN = c95_exp, intercept = FALSE, sf = ama)

munis <- muni_data_bw %>%
    dplyr::filter(inc_study) %>%
    dplyr::arrange(mean) %>%
    dplyr::pull(mun_name) %>%
    unique()

model_labs <- c(
    "(a) Birth-weight mean change",
    "(b) Low birth-weight odds ratio",
    "(c) Prematurity odds ratio"
)
reference <- c(0, 1, 1)

muni_data <- list(muni_data_bw, muni_data_lbw, muni_data_pre) %>%
    map2(model_labs, ~ mutate(.x, model = .y)) %>%
    map2(reference, ~ mutate(.x, reference = .y)) %>%
    bind_rows() %>%
    mutate(mun_name = factor(mun_name, munis))

Credible intervals for municipality random effects

base_size <- 7
ggp_ci <- function(data, munis) {
    data <- data %>%
    dplyr::filter(inc_study) %>%
    dplyr::mutate(mun_name = factor(mun_name, munis))

    gg_ci <- ggplot(data, aes(x = mean, y = mun_name)) +
        geom_errorbarh(aes(xmin = lower, xmax = upper, col = type),
                      position = position_dodge(0.5), height = 0.8, size = rel(0.4)) +
        geom_point(aes(color = type), position = position_dodge(0.5), size = rel(0.7)) +
        geom_vline(aes(xintercept = reference, group = model), lty = 2, col = 1) +
        labs(x = NULL, y = NULL, color = "Controlling gestational age?") +
        facet_wrap(~ model, scale = "free_x") +
        theme_bw(base_size = base_size, base_family = "Helvetica") +
        theme(legend.position = "bottom")
    return(gg_ci)
}

gg_muni_ci <- ggp_ci(muni_data, munis) +
    theme(
          legend.title = element_text(size = base_size, face = "bold"),
          legend.text = element_text(size = base_size),
          axis.text = element_text(size = base_size),
          strip.background = element_blank(),
          strip.text = element_text(size = base_size, hjust = -0.05)
    )

ggsave(file.path(path_summarised, "paper-mod-eff-municipalities-ci.pdf"), gg_muni_ci ,
       width = 18, height = 11.7, units = "cm")
ggsave(file.path(path_summarised, "paper-mod-eff-municipalities-ci.jpg"), gg_muni_ci ,
       width = 18, height = 11.7, units = "cm", dpi = 350)

Visualize municipality random effects with credible intervals

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

Maps of municipality random effects

ggp_map <- function(data, subset, trans = "identity", trans_inv = "identity",
                    rev = FALSE, width = 40, types = c("No", "Yes"), base_size = 7) {

    data <- data[data$model == subset, ]
    mean_max <- max(abs(range(eval(call(trans, data$mean)), na.rm = TRUE)))
    mean_vals <- c(seq(width / 2, mean_max, width), mean_max)
    mean_vals <- eval(call(trans_inv, sort(unique(c(- mean_vals, mean_vals)))))
    data$mean <- cut(data$mean, mean_vals, include.lowest = TRUE)
    mean_levels <- levels(data$mean)
    # data <- within(data, labels(type) <- types)
    # data <- mutate(data)
    levels(data$type) <- types

    pal_args <- list(
        palette = "Purple-Green", l1 = 30, p1 = 0.9, p2 = 1, rev = rev,
        breaks = if (rev) mean_levels else rev(mean_levels),
        na.value = "gray70",
        drop = FALSE
    )

    gg <- data %>%
    ggplot() +
    geom_sf(aes(fill = mean), size = rel(0.2)) +
    do.call(colorspace::scale_fill_discrete_diverging, pal_args) +
    # scale_fill_discrete_diverging(palette = "Purple-Green",
    #                               na.value = "gray70",
    #                               breaks = mean_levels) +
    facet_grid(model ~ type) +
    theme_bw(base_size = base_size, base_family = "Helvetica") +
    theme(strip.background = element_blank(),
          strip.text.x = element_text(hjust = 0, size = base_size),
          strip.text.y = element_blank(),
          legend.key.size = unit(0.7, "lines"),
          plot.margin = margin(0, 0, 0, 0)
          ) +
    labs(fill = NULL)
return(gg)
}

add <- muni_data[1, ]
add$type <- "Yes"
add$model <- model_labs[3]
add$mean <- NA

muni_data2 <- rbind(muni_data, add)

gg_map1 <- ggp_map(muni_data, model_labs[1], types = c("(a)", "(b)")) +
    labs(fill = day2day::wrapit("Birth-weight mean change", 15))
gg_map2 <- ggp_map(muni_data, model_labs[2], "log", "exp",
                   TRUE, 0.1, types = c("(c)", "(d)")) +
    labs(fill = "LBW odds ratio")

blank <- plot_spacer() +
    theme(axis.text.y = element_blank(), plot.margin = margin(0, 0, 0, 4))
gg_map3 <- ggp_map(muni_data, model_labs[3], "log", "exp",
                   TRUE, 0.3, types = c("(e)", "(f)")) +
    labs(fill = "PTB odds ratio") +
    blank +
    plot_layout(guides = "collect")

base_size <- 7
gg_muni_map <- gg_map1 + gg_map2 + gg_map3 +
    plot_layout(nrow = 3, heights = c(1, 1, 1)) &
    theme(
          legend.title = element_text(size = base_size, face = "bold"),
          legend.text = element_text(size = base_size),
          axis.text = element_text(size = base_size)
    )

ggsave(file.path(path_summarised, "paper-mod-eff-municipalities-map.pdf"), gg_muni_map ,
       width = 18, height = 17.4, units = "cm")
ggsave(file.path(path_summarised, "paper-mod-eff-municipalities-map.jpg"), gg_muni_map ,
       width = 18, height = 17.4, units = "cm", dpi = 350)

Visualize maps of municipality random effects

print(gg_muni_map, 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 
#>  76.464   3.017  79.496