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.construct
library(gamlss.dist)
#> Loading required package: MASS
path_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.708

Parametric 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