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, "pre-12-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:
#> ---
#> 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) -5.973756 -15.757256 -4.531461 -2.563529 -0.521
#> sexfemale -0.029916 -0.057292 -0.029638 -0.003081 -0.030
#> born_raceindigenous 0.263663 0.220616 0.262950 0.305008 0.263
#> born_raceNA -0.058985 -0.201707 -0.060789 0.082890 -0.057
#> birth_placeanother health center 0.133869 -0.053145 0.129153 0.329914 0.135
#> birth_placehome 0.211896 0.172611 0.211800 0.251473 0.210
#> birth_placeother 0.354375 0.242680 0.356420 0.479038 0.351
#> birth_placeNA -0.109262 -0.240937 -0.111009 0.019001 -0.109
#> marital_statusmarried -0.012972 -0.061497 -0.013082 0.036139 -0.014
#> marital_statuswidow 0.056703 -0.351630 0.058760 0.471346 0.055
#> marital_statusdivorced -0.164206 -0.572626 -0.158932 0.280315 -0.170
#> marital_statusconsensual union -0.051418 -0.085595 -0.051046 -0.020054 -0.054
#> marital_statusNA -0.013284 -0.143167 -0.013493 0.105666 -0.012
#> study_years1 - 3 -0.088592 -0.154175 -0.088509 -0.023549 -0.088
#> study_years4 - 7 -0.084107 -0.145763 -0.085557 -0.021997 -0.084
#> study_years8 - 11 0.001172 -0.061190 0.001339 0.069380 -0.002
#> study_years> 12 0.039171 -0.051407 0.038395 0.136451 0.038
#> study_yearsNA -0.121149 -0.267747 -0.119502 0.040934 -0.124
#> consult_num1 - 3 0.523229 0.463738 0.522894 0.589800 0.522
#> consult_num4 - 6 0.250917 0.190190 0.249243 0.314311 0.251
#> consult_num> 7 -0.529896 -0.595414 -0.529447 -0.460221 -0.529
#> consult_numNA 0.175304 0.013178 0.175608 0.343119 0.175
#> remoteness 3.007561 -2.102575 2.401481 11.168710 -0.065
#> prop_tap_toilet 5.340236 -1.121827 2.435829 24.963426 0.503
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8286 0.2190 0.9402 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 1.358e+02 6.700e+01 1.252e+02 2.606e+02
#> s(age).alpha 6.981e-01 8.152e-03 8.444e-01 1.000e+00
#> s(age).edf 8.036e+00 7.647e+00 8.031e+00 8.460e+00
#> s(wk_ini).tau21 1.506e+02 7.703e+01 1.412e+02 2.932e+02
#> s(wk_ini).alpha 9.691e-01 8.228e-01 9.949e-01 1.000e+00
#> s(wk_ini).edf 8.970e+00 8.946e+00 8.971e+00 8.985e+00
#> s(rivwk_conception).tau21 1.231e+02 6.463e+01 1.152e+02 2.377e+02
#> s(rivwk_conception).alpha 9.821e-01 8.979e-01 9.999e-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.150e+01 2.789e+01 4.012e+01 6.218e+01
#> s(res_muni).alpha 7.818e-01 2.718e-01 8.710e-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 3.626e+01 2.650e+01 3.541e+01 5.009e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 6.793e-01 4.763e-02 7.395e-01 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 6.729e+01 6.670e+01 6.731e+01 6.778e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 3.899e+01 2.822e+01 3.826e+01 5.357e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 7.531e-01 2.018e-01 8.234e-01 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 6.777e+01 6.738e+01 6.778e+01 6.814e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 3.357e+02 2.393e+02 3.306e+02 4.698e+02
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 2.705e-01 3.135e-06 8.225e-02 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 6.828e+01 6.802e+01 6.829e+01 6.850e+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 = 58490.48
#> ---
#> Optimizer summary:
#> -
#> AICc = 152943.9 edf = 288.4333 logLik = -76183.21
#> logPost = -671247.9 nobs = 291479 runtime = 5083.758Parametric 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.923 0.248 7.240