Summarise (not used) preterm birth (PTB) model: gestational age path


In this script we summarise the results of the PTB model. We also include comparisons between different types of profiles. The resulting table can be seen in section Display summary. These results might not be as good as the results from another model with greater number of MCMC iterations.

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_pre <- file.path(path_modelled, "pre-10-full-re-burned.rds")

bwdata_model <- fst::read_fst(bwdata_file)
model_pre <- readRDS(path_model_pre)

data_points_all <- readRDS(file.path(path_processed, "exposure-reference-points.rds"))
data_summary <- tibble::tribble(~description, ~value,
    "number of birth registration (modelled)", as.character(nrow(bwdata_model)),
    "period of study", "2006 - 2017")

# model summary
bamlss_pre_sum <- summary(model_pre)

data_summary <- tibble::tribble(~description, ~value,
        "proportion of prematurity", round(mean(bwdata_model$premature), 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_pre <- predict(model_pre, 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_pre, 1, function (x) c95(exp(x))))

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

data_summary <- tibble::tribble(~description, ~value,
    "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, .)
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 pre
data_extremes_pre <- data_points_all[[3]] %>%
    mutate(description = paste0("odds of extreme events ", label, " (", description, ")"))
# data_extremes_pre[data_extremes_pre$label == "E", "y"] <- 0.4

ext_eff_pre <- predict_effect(model_pre, data_extremes_pre, "pi",
                          "s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk)",
                          FUN = function(x) c95(exp(x)))

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

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

ext_eff_pre_base <- predict_effect_baseline(
    model_pre, data_extremes_pre, "pi",
    "s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk)",
    FUN = c95_exp_prob)

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

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

events_eff_pre <- predict_effect(model_pre, data_events_pre,
                          "pi", "s(neg_mckee_mean_8wk,pos_mckee_mean_8wk)",
                          FUN = function(x) c95(exp(x)))

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

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

events_eff_pre_base <- predict_effect_baseline(
    model_pre, data_events_pre, "pi",
    "s(neg_mckee_mean_8wk,pos_mckee_mean_8wk)",
    FUN = c95_exp_prob
)

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

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

deviations_eff_pre <- predict_effect(model_pre, data_deviations_pre,
                          "pi", "s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk)",
                          FUN = function(x) c95(exp(x)))

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

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

deviations_eff_pre <- predict_effect_baseline(
    model_pre, data_deviations_pre, "pi",
    "s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk)",
    FUN = c95_exp_prob
)

data_summary  <- data_deviations_pre %>%
    dplyr::mutate(value = ci_character_exceed(deviations_eff_pre)) %>%
    dplyr::select(description, value) %>%
    rbind(data_summary, .)
# remoteness municipality level effects (mean and lbw)
data_summary  <-
    effect_summary_prob_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, 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_pre(model_pre, data.frame(c(167, 700)), term = "s(wk_ini)",
                   label_1 = "2006", label_2 = "2016") %>%
    rbind(data_summary, .)

Display summary

cap <- "Summarise preterm birth (PTB) models."
knitr::kable(data_summary, "html", table.attr = "class=\"table\"",
             caption = cap)
Table 1: Summarise preterm birth (PTB) models.
description value
number of birth registration (modelled) 291479
period of study 2006 - 2017
proportion of prematurity 0.091
odds ratio when born to adolescent, indigenous Amerindian, unmarried mothers that received no formal education, antenatal or obstetric healthcare 4.99 (CI: 4.44, 5.55)
odds change (with respect to baseline) when born to adolescent, indigenous Amerindian, unmarried mothers that received no formal education, antenatal or obstetric healthcare 2.33 (CI: 2.19, 2.47)
odds change (with respect to baseline) when born to mature, non-indigenous, married mothers that received formal education (4-7), antenatal or obstetric healthcare 0.47 (CI: 0.42, 0.52)
odds of extreme events A (None) 1.384 (CI: 1.346, 1.423)
odds of extreme events B (Intense and deficient) 5034.317 (CI: 305.277, 21497.17)
odds of extreme events C (Intense) 306.487 (CI: 22.081, 1228.598)
odds of extreme events E (Moderate intense) 0.004 (CI: 0.001, 0.015)
odds of extreme events D (Deficient) 6.045 (CI: 1.396, 16.38)
extreme effects: odds ratio of A - A (None) 1 (CI: 1, 1); P(diff > 0) = 0
extreme effects: odds ratio of B - A (Intense and deficient) 3638.152 (CI: 223.892, 15560.951); P(diff > 0) = 1
extreme effects: odds ratio of C - A (Intense) 221.828 (CI: 15.701, 894.406); P(diff > 0) = 1
extreme effects: odds ratio of E - A (Moderate intense) 0.003 (CI: 0.001, 0.011); P(diff > 0) = 0
extreme effects: odds ratio of D - A (Deficient) 4.377 (CI: 1.01, 11.872); P(diff > 0) = 0.975
odds of non-extreme effects of A (None) 0.99 (CI: 0.9, 1.09)
odds of non-extreme effects of B (Intense and deficient) 1.05 (CI: 0.62, 1.65)
odds of non-extreme effects of C (Intense) 1.76 (CI: 1.12, 2.61)
odds of non-extreme effects of D (Deficient) 0.44 (CI: 0.18, 0.87)
non-extreme effects: odds ratio of A - A (None) 1 (CI: 1, 1); P(diff > 0) = 0
non-extreme effects: odds ratio of B - A (Intense and deficient) 1.062 (CI: 0.623, 1.675); P(diff > 0) = 0.551
non-extreme effects: odds ratio of C - A (Intense) 1.788 (CI: 1.125, 2.668); P(diff > 0) = 0.994
non-extreme effects: odds ratio of D - A (Deficient) 0.448 (CI: 0.176, 0.926); P(diff > 0) = 0.014
odds of seasonal deviation effects of A (None) 0.8 (CI: 0.74, 0.88)
odds of seasonal deviation effects of B (Intense and deficient) 2.24 (CI: 1.34, 3.54)
odds of seasonal deviation effects of C (Intense) 1.77 (CI: 1.23, 2.45)
odds of seasonal deviation effects of D (Deficient) 10.25 (CI: 5.78, 16.75)
seasonal deviation effects: odds ratio of A - A(None) 1 (CI: 1, 1); P(diff > 0) = 0
seasonal deviation effects: odds ratio of B - A(Intense and deficient) 2.784 (CI: 1.657, 4.455); P(diff > 0) = 1
seasonal deviation effects: odds ratio of C - A(Intense) 2.207 (CI: 1.529, 3.058); P(diff > 0) = 1
seasonal deviation effects: odds ratio of D - A(Deficient) 12.786 (CI: 7.026, 21.488); P(diff > 0) = 1
pre: odds ratio between high remoteness (0.75) and low remoteness (0.25) 1.07 (CI: 0.8, 1.39); P(diff > 0) = 0.66
pre: odds change of low remoteness (0.25) 1.03 (CI: 0.9, 1.18); P(diff > 0) = 0.66
pre: odds change of high remoteness (0.75) 1.11 (CI: 0.72, 1.65); P(diff > 0) = 0.66
pre: odds ratio between low sanitation (0.2) and high sanitation (0.45) 0.75 (CI: 0.57, 0.98); P(diff > 0) = 0.02
pre: odds change of high sanitation (0.45) 1.75 (CI: 1.04, 2.79); P(diff > 0) = 0.98
pre: odds change of low sanitation (0.2) 1.27 (CI: 1.02, 1.58); P(diff > 0) = 0.98
pre: odds ratio between remote and low sanitation (0.8, 0.15) and not remote and high sanitation (0.25, 0.5) 0.72 (CI: 0.5, 1.01); P(diff > 0) = 0.03
pre: odds change of not remote and high sanitation (0.25, 0.5) 1.94 (CI: 1.03, 3.56); P(diff > 0) = 0.98
pre: odds change of remote and low sanitation (0.8, 0.15) 1.35 (CI: 0.78, 2.29); P(diff > 0) = 0.83
pre: odds ratio between young mothers (15) and mature mothers (30) 1.66 (CI: 1.58, 1.74)
pre: odds change of mature mothers (30) 0.87 (CI: 0.85, 0.9)
pre: odds change of young mothers (15) 1.45 (CI: 1.4, 1.49)
pre: odds ratio between mothers with no formal education (no) and mother with education (1 - 3) 1.09 (CI: 1.02, 1.17)
pre: odds change of mother with education (1 - 3) 0.92 (CI: 0.85, 0.98)
pre: odds change of mothers with no formal education (no) 1 (CI: 1, 1)
pre: odds ratio between mothers with no formal education (no) and mother with education (4 - 7) 1.09 (CI: 1.02, 1.16)
pre: odds change of mother with education (4 - 7) 0.92 (CI: 0.86, 0.98)
pre: odds change of mothers with no formal education (no) 1 (CI: 1, 1)
pre: odds ratio between single mothers (single) and married mothers (married) 1.01 (CI: 0.97, 1.07)
pre: odds change of married mothers (married) 0.99 (CI: 0.94, 1.04)
pre: odds change of single mothers (single) 1 (CI: 1, 1)
pre: odds ratio between mothers with no ante-natal care (0) and mothers with ante-natal care (> 7) 1.7 (CI: 1.58, 1.82)
pre: odds change of mothers with ante-natal care (> 7) 0.59 (CI: 0.55, 0.63)
pre: odds change of mothers with no ante-natal care (0) 1 (CI: 1, 1)
pre: odds ratio between births outside medical facility (home) (home) and births in medical facility (hospital) (hospital) 1.23 (CI: 1.19, 1.28)
pre: odds change of births in medical facility (hospital) (hospital) 1 (CI: 1, 1)
pre: odds change of births outside medical facility (home) (home) 1.23 (CI: 1.19, 1.28)
pre: odds ratio between indigenous (indigenous) and non-indigenous (non-indigenous) 1.3 (CI: 1.25, 1.36)
pre: odds change of non-indigenous (non-indigenous) 1 (CI: 1, 1)
pre: odds change of indigenous (indigenous) 1.3 (CI: 1.25, 1.36)
pre: odds ratio between dry season (30) and wet season (0) 0.865 (CI: 0.818, 0.914); P(diff > 0) = 0
pre: odds change of wet season (0) 1.11 (CI: 1.072, 1.148); P(diff > 0) = 1
pre: odds change of dry season (30) 0.96 (CI: 0.924, 0.996); P(diff > 0) = 0.013
pre: odds ratio between rising season (39) and wet season (5) 0.709 (CI: 0.673, 0.746); P(diff > 0) = 0
pre: odds change of wet season (5) 1.168 (CI: 1.131, 1.208); P(diff > 0) = 1
pre: odds change of rising season (39) 0.828 (CI: 0.801, 0.858); P(diff > 0) = 0
pre: odds ratio between 2016 (700) and 2006 (167) 7.44 (CI: 6.1, 8.99)
pre: odds change of 2006 (167) 0.22 (CI: 0.18, 0.26)
pre: odds change of 2016 (700) 1.6 (CI: 1.53, 1.67)

Write summary of the data and model

readr::write_csv(as.data.frame(data_summary),
                 file.path(path_summarised, "summary-prematurity-model-not-used.csv"))

Time to execute the task

Only useful when executed with Rscript.

proc.time()
#>    user  system elapsed 
#>  19.516   1.093  20.798