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.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-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.335

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 
#>   6.942   0.170   7.135