Full model with flexible prior: 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, "lbw-11-full-re-p1.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 = 8000,
#> burnin = 0)
#> ---
#> Family: binomial
#> Link function: pi = logit
#> *---
#> Formula pi:
#> ---
#> lbw ~ sex + born_race + birth_place + marital_status + study_years +
#> consult_num + s(age, xt = hyper) + s(wk_ini, xt = hyper) +
#> s(rivwk_conception, bs = "cc", xt = hyper) + remoteness +
#> prop_tap_toilet + s(res_muni, bs = "re", xt = hyper) + s(neg_mbsi_mean_1wk,
#> pos_mbsi_mean_1wk, k = 70, xt = hyper) + s(neg_mckee_mean_8wk,
#> pos_mckee_mean_8wk, k = 70, xt = hyper) + s(neg_ext_mbsi_mean_8wk,
#> pos_ext_mbsi_mean_8wk, k = 70, xt = hyper)
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) -5.746139 -9.648337 -5.751504 -2.066154 0.593
#> sexfemale 0.206497 0.175750 0.206411 0.238341 0.207
#> born_raceindigenous 0.064123 0.011987 0.064634 0.113365 0.065
#> born_raceNA -0.011340 -0.202184 -0.013212 0.187688 -0.013
#> birth_placeanother health center -0.003873 -0.251103 -0.007934 0.261514 0.000
#> birth_placehome 0.014958 -0.032385 0.014655 0.063820 0.014
#> birth_placeother 0.340737 0.200700 0.341058 0.484380 0.340
#> birth_placeNA 0.024059 -0.140160 0.022798 0.193054 0.040
#> marital_statusmarried -0.159074 -0.215010 -0.158473 -0.106191 -0.160
#> marital_statuswidow -0.326396 -0.868855 -0.329963 0.252033 -0.340
#> marital_statusdivorced -0.126463 -0.587925 -0.128063 0.342259 -0.131
#> marital_statusconsensual union -0.162087 -0.204846 -0.162393 -0.119899 -0.164
#> marital_statusNA 0.027757 -0.105924 0.029033 0.164145 0.026
#> study_years1 - 3 -0.247644 -0.317036 -0.248742 -0.177829 -0.248
#> study_years4 - 7 -0.290916 -0.359794 -0.290158 -0.224740 -0.290
#> study_years8 - 11 -0.206549 -0.275229 -0.205026 -0.139143 -0.208
#> study_years> 12 -0.038393 -0.132483 -0.038768 0.059772 -0.038
#> study_yearsNA -0.228462 -0.419522 -0.228442 -0.035374 -0.233
#> consult_num1 - 3 -0.299061 -0.356183 -0.298951 -0.242306 -0.298
#> consult_num4 - 6 -0.603640 -0.660887 -0.603703 -0.546783 -0.603
#> consult_num> 7 -1.061549 -1.129192 -1.061941 -0.996885 -1.059
#> consult_numNA -0.005738 -0.142005 -0.005864 0.142567 -0.005
#> remoteness 2.756496 0.363372 2.688799 5.171101 -0.041
#> prop_tap_toilet 7.697790 -2.432591 7.975479 17.659263 0.127
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.74508 0.07741 0.84548 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 1.390e+02 7.038e+01 1.280e+02 2.728e+02
#> s(age).alpha 6.547e-01 1.148e-03 7.623e-01 1.000e+00
#> s(age).edf 7.976e+00 7.557e+00 7.962e+00 8.429e+00
#> s(wk_ini).tau21 1.279e+02 6.583e+01 1.176e+02 2.459e+02
#> s(wk_ini).alpha 9.662e-01 8.035e-01 9.968e-01 1.000e+00
#> s(wk_ini).edf 8.965e+00 8.939e+00 8.966e+00 8.983e+00
#> s(rivwk_conception).tau21 1.270e+02 6.377e+01 1.176e+02 2.463e+02
#> s(rivwk_conception).alpha 9.733e-01 8.499e-01 9.980e-01 1.000e+00
#> s(rivwk_conception).edf 8.000e+00 8.000e+00 8.000e+00 8.000e+00
#> s(res_muni).tau21 4.088e+01 2.759e+01 3.995e+01 6.011e+01
#> s(res_muni).alpha 7.229e-01 1.423e-01 7.865e-01 1.000e+00
#> s(res_muni).edf 4.300e+01 4.299e+01 4.300e+01 4.300e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 2.921e+01 2.118e+01 2.881e+01 4.025e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 6.022e-01 2.697e-02 5.993e-01 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 6.660e+01 6.587e+01 6.662e+01 6.727e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 2.962e+01 2.136e+01 2.918e+01 4.086e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 6.578e-01 6.696e-02 6.809e-01 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 6.712e+01 6.651e+01 6.713e+01 6.765e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 3.861e+01 2.750e+01 3.789e+01 5.544e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 3.259e-01 1.176e-04 1.418e-01 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 6.527e+01 6.410e+01 6.529e+01 6.632e+01
#> parameters
#> s(age).tau21 0.690
#> s(age).alpha NA
#> s(age).edf 4.630
#> s(wk_ini).tau21 6.608
#> s(wk_ini).alpha NA
#> s(wk_ini).edf 8.483
#> s(rivwk_conception).tau21 0.000
#> s(rivwk_conception).alpha NA
#> s(rivwk_conception).edf 0.255
#> s(res_muni).tau21 7.536
#> s(res_muni).alpha NA
#> s(res_muni).edf 42.975
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 0.320
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha NA
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 28.628
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 1.538
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha NA
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 49.619
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 21.716
#> 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 63.174
#> ---
#> Sampler summary:
#> -
#> runtime = 27381.24
#> ---
#> Optimizer summary:
#> -
#> AICc = 129526.1 edf = 221.7651 logLik = -64541.13
#> logPost = -1137510546 nobs = 291479 runtime = 610.701Parametric 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, image = TRUE, contour = TRUE)
Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 7.028 0.269 7.365