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.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, "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.758

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, 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