Full model: 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-10-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 = 3000,
#> burnin = 0)
#> ---
#> Family: binomial
#> Link function: pi = logit
#> *---
#> Formula pi:
#> ---
#> lbw ~ 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) -2.1658376 -2.4319624 -2.1601403 -1.9060922 0.593
#> sexfemale 0.2070039 0.1759772 0.2066543 0.2363417 0.207
#> born_raceindigenous 0.0667911 0.0174085 0.0674261 0.1169290 0.065
#> born_raceNA -0.0176114 -0.2073267 -0.0174042 0.1798509 -0.013
#> birth_placeanother health center 0.0011679 -0.2391878 0.0003536 0.2500118 0.000
#> birth_placehome 0.0120588 -0.0354034 0.0116492 0.0594698 0.014
#> birth_placeother 0.3358386 0.1957715 0.3374128 0.4839088 0.340
#> birth_placeNA 0.0528583 -0.1084978 0.0494021 0.2087015 0.040
#> marital_statusmarried -0.1577649 -0.2132233 -0.1581320 -0.1018063 -0.160
#> marital_statuswidow -0.3187990 -0.8401923 -0.3281749 0.2644606 -0.340
#> marital_statusdivorced -0.1238819 -0.5988271 -0.1402106 0.3874013 -0.131
#> marital_statusconsensual union -0.1589283 -0.1987011 -0.1593350 -0.1195457 -0.164
#> marital_statusNA 0.0236716 -0.1203158 0.0254596 0.1624749 0.026
#> study_years1 - 3 -0.2469705 -0.3203908 -0.2466745 -0.1767881 -0.248
#> study_years4 - 7 -0.2878626 -0.3581063 -0.2887641 -0.2204480 -0.290
#> study_years8 - 11 -0.2034260 -0.2739073 -0.2038168 -0.1328680 -0.208
#> study_years> 12 -0.0322366 -0.1322795 -0.0315433 0.0645347 -0.038
#> study_yearsNA -0.2489579 -0.4412563 -0.2483397 -0.0526108 -0.233
#> consult_num1 - 3 -0.2961960 -0.3501537 -0.2961071 -0.2392788 -0.298
#> consult_num4 - 6 -0.6023635 -0.6591067 -0.6029809 -0.5466483 -0.603
#> consult_num> 7 -1.0586371 -1.1235168 -1.0597880 -0.9929174 -1.059
#> consult_numNA -0.0024748 -0.1356399 -0.0032571 0.1392619 -0.005
#> remoteness 0.1968927 -0.0592087 0.1937581 0.4669553 -0.041
#> prop_tap_toilet -0.0111101 -0.4971463 -0.0146280 0.4850063 0.127
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.74246 0.04904 0.84867 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 3.890e+00 1.088e+00 2.806e+00 1.270e+01
#> s(age).alpha 9.008e-01 2.871e-01 9.886e-01 1.000e+00
#> s(age).edf 5.904e+00 5.016e+00 5.865e+00 6.937e+00
#> s(wk_ini).tau21 1.124e+00 1.620e-02 6.996e-01 4.523e+00
#> s(wk_ini).alpha 9.777e-01 8.625e-01 9.991e-01 1.000e+00
#> s(wk_ini).edf 6.469e+00 2.950e+00 6.694e+00 8.340e+00
#> s(rivwk_conception).tau21 1.385e-04 2.552e-05 9.461e-05 5.202e-04
#> s(rivwk_conception).alpha 9.923e-01 9.472e-01 1.000e+00 1.000e+00
#> s(rivwk_conception).edf 3.293e+00 2.147e+00 3.198e+00 5.069e+00
#> s(res_muni).tau21 2.383e-02 1.390e-02 2.284e-02 3.782e-02
#> s(res_muni).alpha 7.901e-01 2.818e-01 8.718e-01 1.000e+00
#> s(res_muni).edf 3.639e+01 3.349e+01 3.648e+01 3.875e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 7.075e-02 3.858e-03 5.049e-02 2.667e-01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 9.310e-01 5.978e-01 9.908e-01 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 1.411e+01 4.415e+00 1.358e+01 2.691e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 2.650e-01 6.851e-02 2.228e-01 6.797e-01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 8.693e-01 4.358e-01 9.496e-01 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 2.812e+01 1.751e+01 2.792e+01 3.946e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 6.021e+00 3.179e+00 5.831e+00 1.017e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 4.501e-01 3.423e-03 3.292e-01 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 5.402e+01 4.859e+01 5.418e+01 5.838e+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 = 10561.22
#> ---
#> Optimizer summary:
#> -
#> AICc = 129526.1 edf = 221.7651 logLik = -64541.13
#> logPost = -65039.6 nobs = 291479 runtime = 606.335Parametric 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
#> 6.942 0.170 7.135