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, "pre-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)Compute results
model$results <- results.bamlss.default(model)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:
#> ---
#> premature ~ 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) -3.391953 -3.977223 -3.387707 -2.902758 -0.521
#> sexfemale -0.030518 -0.056955 -0.030234 -0.004291 -0.030
#> born_raceindigenous 0.264314 0.222301 0.264833 0.306114 0.263
#> born_raceNA -0.055330 -0.203424 -0.055710 0.101515 -0.057
#> birth_placeanother health center 0.134670 -0.062372 0.131683 0.337084 0.135
#> birth_placehome 0.210549 0.170978 0.210427 0.249573 0.210
#> birth_placeother 0.347349 0.231577 0.347833 0.457127 0.351
#> birth_placeNA -0.109367 -0.243168 -0.109187 0.028280 -0.109
#> marital_statusmarried -0.013319 -0.065147 -0.012768 0.034778 -0.014
#> marital_statuswidow 0.068821 -0.347439 0.072311 0.491392 0.055
#> marital_statusdivorced -0.166958 -0.588298 -0.173153 0.286545 -0.170
#> marital_statusconsensual union -0.052523 -0.085329 -0.052435 -0.020035 -0.054
#> marital_statusNA -0.009328 -0.134015 -0.011109 0.115272 -0.012
#> study_years1 - 3 -0.087639 -0.156711 -0.086825 -0.018604 -0.088
#> study_years4 - 7 -0.082565 -0.149097 -0.081279 -0.019123 -0.084
#> study_years8 - 11 0.001068 -0.066717 0.001000 0.069116 -0.002
#> study_years> 12 0.041373 -0.055784 0.042512 0.138367 0.038
#> study_yearsNA -0.121833 -0.270308 -0.121695 0.029983 -0.124
#> consult_num1 - 3 0.523519 0.459660 0.522779 0.589972 0.522
#> consult_num4 - 6 0.252352 0.189157 0.251945 0.317545 0.251
#> consult_num> 7 -0.528730 -0.596559 -0.528847 -0.460115 -0.529
#> consult_numNA 0.176870 0.010936 0.179007 0.339328 0.175
#> remoteness 0.109661 -0.437420 0.111729 0.663971 -0.065
#> prop_tap_toilet 1.175727 0.087228 1.187429 2.283509 0.503
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8257 0.2187 0.9331 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 2.140e+00 4.498e-01 1.487e+00 7.642e+00
#> s(age).alpha 9.310e-01 5.081e-01 9.921e-01 1.000e+00
#> s(age).edf 5.511e+00 4.480e+00 5.503e+00 6.681e+00
#> s(wk_ini).tau21 6.062e+01 2.057e+01 4.961e+01 1.699e+02
#> s(wk_ini).alpha 9.702e-01 8.198e-01 9.979e-01 1.000e+00
#> s(wk_ini).edf 8.918e+00 8.826e+00 8.923e+00 8.979e+00
#> s(rivwk_conception).tau21 2.493e-03 6.634e-04 1.972e-03 7.408e-03
#> s(rivwk_conception).alpha 9.857e-01 9.136e-01 1.000e+00 1.000e+00
#> s(rivwk_conception).edf 6.646e+00 5.619e+00 6.673e+00 7.524e+00
#> s(res_muni).tau21 1.357e-01 8.557e-02 1.307e-01 2.120e-01
#> s(res_muni).alpha 8.090e-01 2.872e-01 9.088e-01 1.000e+00
#> s(res_muni).edf 4.189e+01 4.139e+01 4.191e+01 4.232e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 7.771e+00 5.135e+00 7.517e+00 1.163e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 6.738e-01 4.441e-04 7.462e-01 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 6.230e+01 5.981e+01 6.231e+01 6.435e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 1.175e+01 7.951e+00 1.152e+01 1.762e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 7.605e-01 2.106e-01 8.344e-01 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 6.522e+01 6.383e+01 6.527e+01 6.642e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 3.578e+02 2.518e+02 3.518e+02 5.099e+02
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 2.856e-01 2.220e-06 1.028e-01 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 6.832e+01 6.808e+01 6.833e+01 6.854e+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 = 21055.63
#> ---
#> Optimizer summary:
#> -
#> AICc = 152943.9 edf = 288.4333 logLik = -76183.21
#> logPost = -76982.72 nobs = 291479 runtime = 5075.584Parametric 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
#> 233.840 1.872 235.813