Municipality random 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(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