Summarise BW and LBW models without controlling for gestational age: growth or gestational age path

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

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", "50-summarization.R"))
source(file.path(path_proj, "src", "51-bamlss.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_lbw_full <- file.path(path_modelled, "lbw-10-full-re-burned.rds")

bwdata_model <- fst::read_fst(bwdata_file)
model_bw <- readRDS(path_model_bw_full)
model_lbw <- readRDS(path_model_lbw_full)

data_points_all <- readRDS(file.path(path_processed, "exposure-reference-points.rds"))
ama <- sf::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

data_summary <- tibble::tribble(~description, ~value,
    "number of birth registration (modelled)", as.character(nrow(bwdata_model)),
    "period of study", "2006 - 2017")

# model summary
bamlss_sum <- summary(model_bw)
bamlss_lbw_sum <- summary(model_lbw)

data_summary <- tibble::tribble(~description, ~value,
        "mean birthweight", round(mean(bwdata_model$born_weight), 3),
        "mean birthweight on baseline group", round(bamlss_sum$model.matrix$mu[1, 1], 3),
        "proportion of low birthweight", round(mean(bwdata_model$born_weight < 2500), 3)
        ) %>%
    rbind(data_summary, .)

# disadvantage profiles
profiles <- data.frame(
  age = c(15, 30),
  born_race = factor(c("indigenous", "non-indigenous"),
                     levels = levels(bwdata_model$born_race)),
  study_years = factor(c("no", "4 - 7"),
                          levels = levels(bwdata_model$study_years)),
  marital_status = factor(c("single", "married"),
                          levels = levels(bwdata_model$marital_status)),
  consult_num = factor(c("0", "> 7"),
                          levels = levels(bwdata_model$consult_num)),
  birth_place = factor(c("home", "hospital"),
                          levels = levels(bwdata_model$birth_place))
  )

disadvantage_eff <- predict(model_bw, profiles, "mu",
        c("s(age)", "born_race", "study_years", "marital_status", "consult_num", "birth_place"),
        intercept = FALSE , FUN = function (x) x) %>% as.matrix() %>% diff() %>% as.numeric() %>% c95()

disadvantage_eff_lbw <- predict(model_lbw, profiles, "pi",
        c("s(age)", "born_race", "study_years", "marital_status", "consult_num", "birth_place"),
        intercept = FALSE, FUN = function (x) x)

odds <- t(apply(disadvantage_eff_lbw, 1, function (x) c95(exp(x))))

disadvantage_odds_ratio <- disadvantage_eff_lbw[2:1, ] %>% as.matrix() %>% diff() %>% as.numeric() %>% exp() %>% c95()

data_summary <- tibble::tribble(~description, ~value,
    "penalty when born to adolescent, indigenous Amerindian, unmarried mothers that received no formal education, antenatal or obstetric healthcare", ci_character(disadvantage_eff),
    "odds ratio when born to adolescent, indigenous Amerindian, unmarried mothers that received no formal education, antenatal or obstetric healthcare", ci_character(disadvantage_odds_ratio),
    "odds change (with respect to baseline) when born to adolescent, indigenous Amerindian, unmarried mothers that received no formal education, antenatal or obstetric healthcare", ci_character(odds[1, ]),
    "odds change (with respect to baseline) when born to mature, non-indigenous, married mothers that received formal education (4-7), antenatal or obstetric healthcare", ci_character(odds[2, ])
    ) %>%
    rbind(data_summary, .)

profiles_advantage <- data.frame(
  age = rep(30, 2),
  sex = factor(c("male", "female"),
                     levels = levels(bwdata_model$sex)),
  born_race = factor(rep("non-indigenous", 2),
                     levels = levels(bwdata_model$born_race)),
  study_years = factor(rep("4 - 7", 2),
                          levels = levels(bwdata_model$study_years)),
  marital_status = factor(rep("married", 2),
                          levels = levels(bwdata_model$marital_status)),
  consult_num = factor(rep("> 7", 2),
                          levels = levels(bwdata_model$consult_num)),
  birth_place = factor(rep("hospital", 2),
                          levels = levels(bwdata_model$birth_place))
  )

profiles_advantage_mean <- predict(model_bw, profiles_advantage, "mu",
        c("s(age)", "sex", "born_race", "study_years", "marital_status", "consult_num", "birth_place"),
        intercept = TRUE , FUN = function (x) c95(x))

c95_prob <- function(x)  c(c95(x), exceed = mean(x > 0))
c95_exp_prob <- function(x) c(c95(exp(x)), exceed = mean(x > 0))

# extremes bw
data_extremes <- data_points_all[[3]] %>%
    mutate(description = paste0(" extreme effects of ", label, " (", description, ")"))

ext_eff <- predict_effect(model_bw, data_extremes,
                          "mu", "s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk)")

data_summary  <-
    data_extremes %>%
    dplyr::mutate(value =   ci_character(ext_eff)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

data_extremes <- data_points_all[[3]] %>%
    mutate(description = paste0("extreme effects: baseline comparison of ",
                                label, " - A (", description, ")"))

ext_eff_base <- predict_effect_baseline(
    model_bw, data_extremes,
    "mu", "s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk)", FUN = c95_prob
)

data_summary  <-
    data_extremes %>%
    dplyr::mutate(value =   ci_character_exceed(ext_eff_base)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

# extremes lbw
data_extremes_lbw <- data_points_all[[3]] %>%
    mutate(description = paste0("odds of extreme events ", label, " (", description, ")"))
data_extremes_lbw[data_extremes_lbw$label == "E", "y"] <- 0.4

ext_eff_lbw <- predict_effect(model_lbw, data_extremes_lbw, "pi",
                          "s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk)",
                          FUN = function(x) c95(exp(x)))

data_summary  <- data_extremes_lbw %>%
    dplyr::mutate(value = ci_character(ext_eff_lbw, 3)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

data_extremes_lbw <- data_points_all[[3]] %>%
    mutate(description = paste0("extreme effects: odds ratio of ",
                                label, " - A (", description, ")"))
data_extremes_lbw[data_extremes_lbw$label == "E", "y"] <- 0.4

ext_eff_lbw_base <- predict_effect_baseline(
    model_lbw, data_extremes_lbw, "pi",
    "s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk)",
    FUN = c95_exp_prob)

data_summary  <- data_extremes_lbw %>%
    dplyr::mutate(value = ci_character_exceed(ext_eff_lbw_base, 3)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

# non extreme events
data_events <- data_points_all[[2]] %>%
    mutate(description = paste0("non-extreme effects of ", label, " (", description, ")"))

events_eff <- predict_effect(model_bw, data_events,
                          "mu", "s(neg_mckee_mean_8wk,pos_mckee_mean_8wk)")

data_summary  <- data_events %>%
    dplyr::mutate(value = ci_character(events_eff)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

data_events <- data_points_all[[2]] %>%
    mutate(description = paste0("non-extreme effects: baseline comparison of ",
                                label, " - A (", description, ")"))

events_eff_base <- predict_effect_baseline(
    model_bw, data_events, "mu",
    "s(neg_mckee_mean_8wk,pos_mckee_mean_8wk)",
    FUN = c95_prob
)

data_summary  <- data_events %>%
    dplyr::mutate(value = ci_character_exceed(events_eff_base)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

# non extreme lbw
data_events_lbw <- data_points_all[[2]] %>%
    mutate(description = paste0("odds of non-extreme effects of ",
                                label, " (", description, ")"))

events_eff_lbw <- predict_effect(model_lbw, data_events_lbw,
                          "pi", "s(neg_mckee_mean_8wk,pos_mckee_mean_8wk)",
                          FUN = function(x) c95(exp(x)))

data_summary  <- data_events_lbw %>%
    dplyr::mutate(value = ci_character(events_eff_lbw)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

data_events_lbw <- data_points_all[[2]] %>%
    mutate(description = paste0("non-extreme effects: odds ratio of ",
                                label, " - A (", description, ")"))

events_eff_lbw_base <- predict_effect_baseline(
    model_lbw, data_events_lbw, "pi",
    "s(neg_mckee_mean_8wk,pos_mckee_mean_8wk)",
    FUN = c95_exp_prob)

data_summary  <- data_events_lbw %>%
    dplyr::mutate(value = ci_character_exceed(events_eff_lbw_base)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

# deviations from seasonality
data_deviations <- data_points_all[[1]] %>%
    mutate(description = paste0("seasonal deviation effects of ", label, " (",
                                description, ")"))

deviations_eff <- predict_effect(model_bw, data_deviations,
                          "mu", "s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk)")

data_summary  <- data_deviations %>%
    dplyr::mutate(value = ci_character(deviations_eff)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

data_deviations <- data_points_all[[1]] %>%
    mutate(description = paste0("seasonal deviation effects: baseline comparison of ",
                                label, " - A (", description, ")"))

deviations_eff <- predict_effect_baseline(
    model_bw, data_deviations, "mu",
    "s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk)",
    FUN = c95_prob
)

data_summary  <- data_deviations %>%
    dplyr::mutate(value = ci_character_exceed(deviations_eff)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

# deviations from seasonality lbw
data_deviations_lbw <- data_points_all[[1]] %>%
    mutate(description = paste0("odds of seasonal deviation effects of ", label, " (",
                                description, ")"))

deviations_eff_lbw <- predict_effect(model_lbw, data_deviations_lbw,
                          "pi", "s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk)",
                          FUN = function(x) c95(exp(x)))

data_summary  <- data_deviations_lbw %>%
    dplyr::mutate(value = ci_character(deviations_eff_lbw)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

data_deviations_lbw <- data_points_all[[1]] %>%
    mutate(description = paste0("seasonal deviation effects: odds ratio of ",
                                label, " - A(", description, ")"))

deviations_eff_lbw <- predict_effect_baseline(
    model_lbw, data_deviations_lbw, "pi",
    "s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk)",
    FUN = c95_exp_prob)

data_summary  <- data_deviations_lbw %>%
    dplyr::mutate(value = ci_character_exceed(deviations_eff_lbw)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)

# remoteness municipality level effects (mean and lbw)
data_summary  <-
    effect_summary_prob(model_bw, model_lbw, data.frame(c(0.25, 0.75)),
                   term = "remoteness",
                   label_1 = "low remoteness", label_2 = "high remoteness") %>%
    rbind(data_summary, .)

# sanitation municipality level effects (mean and lbw)
data_summary  <-
    effect_summary_prob(model_bw, model_lbw, data.frame(c(0.45, 0.20)),
                   term = "prop_tap_toilet",
                   label_1 = "high sanitation", label_2 = "low sanitation") %>%
    rbind(data_summary, .)

# municipality profiles
profiles_mun <- tibble::tribble(~x, ~y,
                                0.25, 0.50,
                                0.8, 0.15)
data_summary  <-
    effect_summary_prob(model_bw, model_lbw, profiles_mun,
                   term = c("remoteness",  "prop_tap_toilet"),
                   label_1 = "not remote and high sanitation",
                   label_2 = "remote and low sanitation") %>%
    rbind(data_summary, .)

# age effect (mean and lbw)
data_summary  <-
    effect_summary(model_bw, model_lbw, data.frame(c(30, 15)), term = "s(age)",
                   label_1 = "mature mothers", label_2 = "young mothers") %>%
    rbind(data_summary, .)

# education effect (mean and lbw)
data_study <- data.frame(factor(c("1 - 3", "no"), levels(bwdata_model$study_years)))
data_summary  <-
    effect_summary(model_bw, model_lbw, data_study, term = "study_years",
                   label_1 = "mother with education", label_2 = "mothers with no formal education") %>%
    rbind(data_summary, .)

# education effect (mean and lbw)
data_aux <- data.frame(factor(c("4 - 7", "no"), levels(bwdata_model$study_years)))
data_summary  <-
    effect_summary(model_bw, model_lbw, data_aux, term = "study_years",
                   label_1 = "mother with education", label_2 = "mothers with no formal education") %>%
    rbind(data_summary, .)

# marital status effect (mean and lbw)
data_aux <- data.frame(factor(c("married", "single"), levels(bwdata_model$marital_status)))
data_summary  <-
    effect_summary(model_bw, model_lbw, data_aux, term = "marital_status",
                   label_1 = "married mothers", label_2 = "single mothers") %>%
    rbind(data_summary, .)

# antenatal care effect (mean and lbw)
data_aux <- data.frame(factor(c("> 7", "0"), levels(bwdata_model$consult_num)))
data_summary  <-
    effect_summary(model_bw, model_lbw, data_aux, term = "consult_num",
                   label_1 = "mothers with ante-natal care", label_2 = "mothers with no ante-natal care") %>%
    rbind(data_summary, .)

# birth place effect (mean and lbw)
data_aux <- data.frame(factor(c("hospital", "home"), levels(bwdata_model$birth_place)))
data_summary  <-
    effect_summary(model_bw, model_lbw, data_aux, term = "birth_place",
                   label_1 = "births in medical facility (hospital)", label_2 = "births outside medical facility (home)") %>%
    rbind(data_summary, .)

# indigenous effect (mean and lbw)
data_aux <- data.frame(factor(c("non-indigenous", "indigenous"), levels(bwdata_model$born_race)))
data_summary  <-
    effect_summary(model_bw, model_lbw, data_aux, term = "born_race",
                   label_1 = "non-indigenous", label_2 = "indigenous") %>%
    rbind(data_summary, .)

# river week seasonal effect (mean and lbw)
data_summary  <-
    effect_summary_prob(model_bw, model_lbw, data.frame(c(0, 30)), term = "s(rivwk_conception)",
                   label_1 = "wet season", label_2 = "dry season", digits = 3) %>%
    rbind(data_summary, .)

# river week seasonal effect (mean and lbw)
data_summary  <-
    effect_summary_prob(model_bw, model_lbw, data.frame(c(5, 39)), term = "s(rivwk_conception)",
                   label_1 = "wet season", label_2 = "rising season", digits = 3) %>%
    rbind(data_summary, .)


# temporal effect (mean and lbw)
data_summary  <-
    effect_summary(model_bw, model_lbw, data.frame(c(167, 700)), term = "s(wk_ini)",
                   label_1 = "2006", label_2 = "2016") %>%
    rbind(data_summary, .)

Municipality comparisons

mun_names <- c("Envira", "Carauari")
mun_codes <- ama$mun_code[match(mun_names, ama$mun_name)]
data_aux <- data.frame(factor(mun_codes, levels(bwdata_model$res_muni)))

data_summary  <-
    effect_summary(model_bw, model_lbw, data_aux, term = "s(res_muni)",
               label_1 = mun_names[1], label_2 = mun_names[2]) %>%
    rbind(data_summary, .)

Display summary

cap <- "Summarise BW and LBW models without controlling for gestational age."
knitr::kable(data_summary, "html", table.attr = "class=\"table\"",
             caption = cap)

Write summary of the data and model

readr::write_csv(as.data.frame(data_summary),
                 file.path(path_summarised, "summary-full-models.csv"))

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  32.868   1.640  34.485