Full model with smoother 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-12-full-re-p2.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) -2.159511 -2.346907 -2.160038 -1.969596 0.593
#> sexfemale 0.205646 0.173880 0.205608 0.235656 0.207
#> born_raceindigenous 0.068350 0.018692 0.068271 0.117029 0.065
#> born_raceNA -0.013663 -0.202292 -0.014641 0.181295 -0.013
#> birth_placeanother health center 0.005091 -0.236661 0.000719 0.273344 0.000
#> birth_placehome 0.006050 -0.042074 0.006426 0.054249 0.014
#> birth_placeother 0.319448 0.177596 0.318570 0.464291 0.340
#> birth_placeNA 0.046566 -0.108925 0.045370 0.205595 0.040
#> marital_statusmarried -0.158585 -0.213160 -0.157665 -0.104998 -0.160
#> marital_statuswidow -0.329066 -0.873007 -0.339134 0.260659 -0.340
#> marital_statusdivorced -0.140722 -0.589123 -0.148493 0.331590 -0.131
#> marital_statusconsensual union -0.157783 -0.196516 -0.158116 -0.117997 -0.164
#> marital_statusNA 0.012848 -0.128975 0.013552 0.152456 0.026
#> study_years1 - 3 -0.243038 -0.312327 -0.242867 -0.173189 -0.248
#> study_years4 - 7 -0.286147 -0.350998 -0.286374 -0.220371 -0.290
#> study_years8 - 11 -0.206076 -0.271367 -0.205592 -0.138289 -0.208
#> study_years> 12 -0.033720 -0.131942 -0.034731 0.063700 -0.038
#> study_yearsNA -0.270199 -0.460423 -0.272043 -0.080995 -0.233
#> consult_num1 - 3 -0.297799 -0.357649 -0.298008 -0.239170 -0.298
#> consult_num4 - 6 -0.606966 -0.665467 -0.607325 -0.548905 -0.603
#> consult_num> 7 -1.065997 -1.127078 -1.066704 -1.000670 -1.059
#> consult_numNA -0.009763 -0.144002 -0.011375 0.133401 -0.005
#> remoteness 0.223100 0.022266 0.224474 0.408797 -0.041
#> prop_tap_toilet 0.005369 -0.371073 0.005387 0.382848 0.127
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.73100 0.05576 0.83875 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 5.006e-01 2.731e-01 4.714e-01 9.022e-01
#> s(age).alpha 9.502e-01 6.741e-01 9.936e-01 1.000e+00
#> s(age).edf 4.316e+00 3.806e+00 4.308e+00 4.890e+00
#> s(wk_ini).tau21 1.109e-05 5.921e-06 1.019e-05 2.006e-05
#> s(wk_ini).alpha 9.995e-01 9.965e-01 1.000e+00 1.000e+00
#> s(wk_ini).edf 1.008e+00 1.004e+00 1.007e+00 1.014e+00
#> s(rivwk_conception).tau21 1.075e-05 5.872e-06 1.006e-05 1.925e-05
#> s(rivwk_conception).alpha 9.982e-01 9.875e-01 1.000e+00 1.000e+00
#> s(rivwk_conception).edf 1.443e+00 1.067e+00 1.431e+00 1.885e+00
#> s(res_muni).tau21 1.276e-02 7.856e-03 1.233e-02 1.986e-02
#> s(res_muni).alpha 8.164e-01 3.223e-01 9.043e-01 1.000e+00
#> s(res_muni).edf 3.256e+01 2.906e+01 3.259e+01 3.567e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 1.128e-05 5.963e-06 1.055e-05 2.113e-05
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 9.967e-01 9.755e-01 1.000e+00 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 2.011e+00 2.006e+00 2.010e+00 2.021e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 1.134e-05 5.919e-06 1.051e-05 2.210e-05
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 9.944e-01 9.539e-01 1.000e+00 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 2.017e+00 2.009e+00 2.016e+00 2.032e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 1.102e-05 5.938e-06 1.018e-05 2.094e-05
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 9.840e-01 8.776e-01 9.999e-01 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 2.018e+00 2.010e+00 2.017e+00 2.035e+00
#> 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 = 27525.39
#> ---
#> Optimizer summary:
#> -
#> AICc = 129526.1 edf = 221.7651 logLik = -64541.13
#> logPost = -65629.06 nobs = 291479 runtime = 602.708Parametric 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.030 0.252 7.301