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, "pre-13-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:
#> ---
#> premature ~ 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) -3.321709 -3.710898 -3.327458 -2.890282 -0.521
#> sexfemale -0.030548 -0.058354 -0.030423 -0.004605 -0.030
#> born_raceindigenous 0.265587 0.225745 0.264884 0.310086 0.263
#> born_raceNA -0.049139 -0.191550 -0.051221 0.093767 -0.057
#> birth_placeanother health center 0.145411 -0.060407 0.143541 0.349025 0.135
#> birth_placehome 0.210457 0.172689 0.210390 0.249193 0.210
#> birth_placeother 0.349841 0.232386 0.348696 0.471600 0.351
#> birth_placeNA -0.109602 -0.243538 -0.106660 0.021490 -0.109
#> marital_statusmarried -0.014705 -0.062613 -0.015145 0.035176 -0.014
#> marital_statuswidow 0.040719 -0.376169 0.043265 0.445977 0.055
#> marital_statusdivorced -0.161983 -0.582619 -0.162313 0.301910 -0.170
#> marital_statusconsensual union -0.054799 -0.087593 -0.054575 -0.023186 -0.054
#> marital_statusNA -0.004588 -0.126614 -0.006015 0.125783 -0.012
#> study_years1 - 3 -0.086898 -0.157163 -0.087168 -0.025076 -0.088
#> study_years4 - 7 -0.082733 -0.145568 -0.082317 -0.025400 -0.084
#> study_years8 - 11 -0.002441 -0.062915 -0.001866 0.056117 -0.002
#> study_years> 12 0.037920 -0.061961 0.040030 0.120730 0.038
#> study_yearsNA -0.120912 -0.270881 -0.123257 0.033479 -0.124
#> consult_num1 - 3 0.526183 0.459091 0.526396 0.590791 0.522
#> consult_num4 - 6 0.254413 0.187710 0.254962 0.320474 0.251
#> consult_num> 7 -0.525485 -0.597981 -0.525954 -0.456545 -0.529
#> consult_numNA 0.177374 0.027391 0.177372 0.326808 0.175
#> remoteness 0.051599 -0.420176 0.058109 0.482059 -0.065
#> prop_tap_toilet 1.090935 0.240525 1.104997 1.873832 0.503
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8217 0.2154 0.9197 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 2.017e-01 1.048e-01 1.902e-01 3.799e-01
#> s(age).alpha 9.775e-01 8.339e-01 9.994e-01 1.000e+00
#> s(age).edf 3.699e+00 3.209e+00 3.679e+00 4.270e+00
#> s(wk_ini).tau21 1.349e+01 7.596e+00 1.283e+01 2.434e+01
#> s(wk_ini).alpha 9.714e-01 8.151e-01 9.992e-01 1.000e+00
#> s(wk_ini).edf 8.717e+00 8.568e+00 8.721e+00 8.840e+00
#> s(rivwk_conception).tau21 2.036e-04 7.592e-05 1.854e-04 4.526e-04
#> s(rivwk_conception).alpha 9.918e-01 9.490e-01 9.996e-01 1.000e+00
#> s(rivwk_conception).edf 4.155e+00 3.294e+00 4.155e+00 5.083e+00
#> s(res_muni).tau21 8.614e-02 5.941e-02 8.401e-02 1.241e-01
#> s(res_muni).alpha 8.204e-01 2.956e-01 9.332e-01 1.000e+00
#> s(res_muni).edf 4.133e+01 4.066e+01 4.135e+01 4.187e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 5.067e+00 3.399e+00 4.925e+00 7.414e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 7.274e-01 5.386e-02 8.038e-01 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 5.968e+01 5.677e+01 5.977e+01 6.217e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 8.324e+00 5.697e+00 8.189e+00 1.155e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 7.560e-01 2.243e-01 8.165e-01 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 6.391e+01 6.217e+01 6.397e+01 6.529e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 2.705e+02 2.016e+02 2.690e+02 3.619e+02
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 2.742e-01 2.265e-07 8.118e-02 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 6.813e+01 6.784e+01 6.814e+01 6.835e+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 = 56674.06
#> ---
#> Optimizer summary:
#> -
#> AICc = 152943.9 edf = 288.4333 logLik = -76183.21
#> logPost = -77764.79 nobs = 291479 runtime = 5078.625Parametric 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
#> 6.842 0.299 7.210