Full model 2: effects
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.constructlibrary(gamlss.dist)#> Loading required package: MASSpath_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")
source(file.path(path_proj, "src", "51-bamlss.R"))
bwdata_file <- file.path(path_processed, "bwdata_41_model.fst")
model_file <- file.path(path_modelled, "pre-11-full-re.rds")
form_file <- gsub("(\\.rds)$", "-form\\1", model_file)
model_file_burned <- gsub("(\\.rds)$", "-burned\\1", model_file)
bwdata_model <- fst::read_fst(bwdata_file)
form <- readRDS(form_file)
model <- readRDS(model_file_burned)Summary
summary(model)#>
#> Call:
#> bamlss(formula = form, family = "binomial", data = bwdata_model,
#> cores = 4, combine = FALSE, light = TRUE, n.iter = 10000,
#> burnin = 0)
#> ---
#> Family: binomial
#> Link function: pi = logit
#> *---
#> Formula pi:
#> ---
#> premature ~ sex + born_race + birth_place + marital_status +
#> study_years + consult_num + s(age) + s(wk_ini) + s(rivwk_conception,
#> bs = "cc") + remoteness + prop_tap_toilet + s(res_muni, bs = "re") +
#> s(neg_mbsi_mean_1wk, pos_mbsi_mean_1wk, k = 70) + s(neg_mckee_mean_8wk,
#> pos_mckee_mean_8wk, k = 70) + s(neg_ext_mbsi_mean_8wk, pos_ext_mbsi_mean_8wk,
#> k = 70)
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) -3.346304 -3.851995 -3.342195 -2.878428 -0.521
#> sexfemale -0.030856 -0.056278 -0.030887 -0.005092 -0.030
#> born_raceindigenous 0.264617 0.221416 0.265536 0.304751 0.263
#> born_raceNA -0.053738 -0.196474 -0.054493 0.094917 -0.057
#> birth_placeanother health center 0.142891 -0.054364 0.140934 0.340345 0.135
#> birth_placehome 0.210770 0.170351 0.210940 0.249384 0.210
#> birth_placeother 0.348939 0.227584 0.348386 0.467515 0.351
#> birth_placeNA -0.113087 -0.238591 -0.115099 0.018193 -0.109
#> marital_statusmarried -0.013846 -0.063275 -0.014077 0.039512 -0.014
#> marital_statuswidow 0.055461 -0.371509 0.059530 0.490501 0.055
#> marital_statusdivorced -0.173486 -0.598007 -0.177235 0.283082 -0.170
#> marital_statusconsensual union -0.052287 -0.084219 -0.052148 -0.018822 -0.054
#> marital_statusNA -0.006157 -0.144182 -0.006766 0.121320 -0.012
#> study_years1 - 3 -0.086434 -0.156815 -0.087717 -0.018135 -0.088
#> study_years4 - 7 -0.081894 -0.144564 -0.081489 -0.020704 -0.084
#> study_years8 - 11 0.001122 -0.065013 0.002076 0.062953 -0.002
#> study_years> 12 0.042213 -0.048243 0.041978 0.134227 0.038
#> study_yearsNA -0.121716 -0.267093 -0.123112 0.025241 -0.124
#> consult_num1 - 3 0.523842 0.459627 0.523508 0.587652 0.522
#> consult_num4 - 6 0.251604 0.186033 0.251641 0.318336 0.251
#> consult_num> 7 -0.527994 -0.596431 -0.528202 -0.461309 -0.529
#> consult_numNA 0.172449 0.011358 0.172846 0.330750 0.175
#> remoteness 0.028863 -0.501809 0.029396 0.586696 -0.065
#> prop_tap_toilet 1.099025 0.151872 1.089929 2.107041 0.503
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8140 0.1911 0.9250 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 2.130e+00 4.436e-01 1.460e+00 7.541e+00
#> s(age).alpha 9.359e-01 4.674e-01 9.979e-01 1.000e+00
#> s(age).edf 5.490e+00 4.447e+00 5.454e+00 6.672e+00
#> s(wk_ini).tau21 6.307e+01 2.160e+01 5.211e+01 1.683e+02
#> s(wk_ini).alpha 9.703e-01 7.921e-01 9.990e-01 1.000e+00
#> s(wk_ini).edf 8.922e+00 8.832e+00 8.927e+00 8.978e+00
#> s(rivwk_conception).tau21 2.505e-03 7.120e-04 1.942e-03 7.668e-03
#> s(rivwk_conception).alpha 9.858e-01 9.130e-01 1.000e+00 1.000e+00
#> s(rivwk_conception).edf 6.683e+00 5.737e+00 6.703e+00 7.561e+00
#> s(res_muni).tau21 1.378e-01 8.522e-02 1.316e-01 2.183e-01
#> s(res_muni).alpha 8.041e-01 2.636e-01 8.975e-01 1.000e+00
#> s(res_muni).edf 4.192e+01 4.139e+01 4.193e+01 4.234e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 7.834e+00 5.138e+00 7.599e+00 1.187e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 6.894e-01 1.273e-02 7.632e-01 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 6.236e+01 6.009e+01 6.241e+01 6.446e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 1.183e+01 8.295e+00 1.153e+01 1.744e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 7.594e-01 2.245e-01 8.279e-01 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 6.524e+01 6.399e+01 6.525e+01 6.642e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 3.601e+02 2.527e+02 3.526e+02 4.984e+02
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 2.945e-01 2.105e-05 9.667e-02 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 6.833e+01 6.809e+01 6.834e+01 6.854e+01
#> parameters
#> s(age).tau21 0.424
#> s(age).alpha NA
#> s(age).edf 4.385
#> s(wk_ini).tau21 44.459
#> s(wk_ini).alpha NA
#> s(wk_ini).edf 8.913
#> s(rivwk_conception).tau21 0.002
#> s(rivwk_conception).alpha NA
#> s(rivwk_conception).edf 6.539
#> s(res_muni).tau21 7.365
#> s(res_muni).alpha NA
#> s(res_muni).edf 42.980
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 19.898
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha NA
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 66.104
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 20.700
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha NA
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 66.821
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 772.420
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha NA
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 68.690
#> ---
#> Sampler summary:
#> -
#> runtime = 71308.86
#> ---
#> Optimizer summary:
#> -
#> AICc = 152943.9 edf = 288.4333 logLik = -76183.21
#> logPost = -76982.72 nobs = 291479 runtime = 5078.429Parametric effects
par(mar = c(4, 4, 0.5, 0), mfrow = c(1, 2), cex.axis = 0.7)
plot2d_bamlss(model, bwdata_model, model = "pi", term = "remoteness", grid = 50,
FUN = c95)
plot2d_bamlss(model, bwdata_model, model = "pi", term = "prop_tap_toilet", grid = 50,
FUN = c95)
Smooth effects
par(mar = c(4, 4, 0.5, 0), mfrow = c(1, 2), cex.axis = 0.7)
plot(model, scale = 0, scheme = 2, spar = FALSE, contour = TRUE, image = TRUE)
Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 7.054 0.231 7.304