Gaussian 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, "bw-12-full-re.rds")
model_file_burned <- gsub("(\\.rds)$", "-burned\\1", model_file)
bwdata_model <- fst::read_fst(bwdata_file)
model <- readRDS(model_file_burned)Summary
summary(model)#>
#> Call:
#> bamlss(formula = form, data = bwdata_model, cores = 4, combine = FALSE,
#> light = TRUE, n.iter = 7000, burnin = 0)
#> ---
#> Family: gaussian
#> Link function: mu = identity, sigma = log
#> *---
#> Formula mu:
#> ---
#> born_weight ~ 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) 3077.619 3018.176 3075.890 3139.415 3106.391
#> sexfemale -98.479 -101.988 -98.539 -94.971 -98.517
#> born_raceindigenous -53.840 -59.400 -53.895 -48.534 -53.640
#> born_raceNA 8.850 -12.769 8.526 31.249 8.611
#> birth_placeanother health center 13.077 -18.675 14.086 43.333 12.718
#> birth_placehome -59.856 -65.679 -59.817 -53.710 -59.795
#> birth_placeother -94.457 -115.389 -94.351 -72.897 -94.239
#> birth_placeNA -103.956 -123.952 -104.190 -81.944 -103.600
#> marital_statusmarried 33.531 27.463 33.735 39.415 33.640
#> marital_statuswidow -3.584 -62.977 -3.324 57.171 -4.970
#> marital_statusdivorced -23.704 -74.495 -23.254 28.521 -24.805
#> marital_statusconsensual union 29.965 25.447 29.968 34.798 30.259
#> marital_statusNA 21.871 6.706 21.880 37.624 22.101
#> study_years1 - 3 63.266 53.286 63.386 72.699 63.280
#> study_years4 - 7 81.834 72.990 81.850 90.401 81.951
#> study_years8 - 11 71.959 62.539 71.808 81.233 72.151
#> study_years> 12 65.458 54.005 65.601 77.196 65.876
#> study_yearsNA 47.014 23.689 47.053 69.493 46.507
#> consult_num1 - 3 57.371 49.060 57.295 66.653 57.662
#> consult_num4 - 6 109.385 101.257 109.338 118.302 109.492
#> consult_num> 7 177.518 168.676 177.571 186.472 177.542
#> consult_numNA 7.029 -14.717 7.025 30.339 6.508
#> remoteness -19.046 -85.957 -19.393 47.651 -24.895
#> prop_tap_toilet 98.753 -32.020 101.156 233.388 -90.782
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.9886 0.9446 1.0000 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 1.977e+05 5.542e+04 1.537e+05 5.985e+05
#> s(age).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
#> s(age).edf 6.580e+00 5.855e+00 6.572e+00 7.313e+00
#> s(wk_ini).tau21 3.627e+03 8.268e+01 8.355e+02 3.224e+04
#> s(wk_ini).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
#> s(wk_ini).edf 4.332e+00 2.265e+00 4.031e+00 7.875e+00
#> s(rivwk_conception).tau21 9.478e-01 2.836e-02 4.640e-01 5.182e+00
#> s(rivwk_conception).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
#> s(rivwk_conception).edf 2.390e+00 3.961e-01 2.345e+00 4.623e+00
#> s(res_muni).tau21 1.830e+03 1.167e+03 1.756e+03 2.877e+03
#> s(res_muni).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
#> s(res_muni).edf 4.174e+01 4.115e+01 4.176e+01 4.222e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 1.520e+03 1.041e+02 1.102e+03 4.817e+03
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 1.735e+01 5.853e+00 1.683e+01 3.073e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 3.406e+03 1.183e+03 3.041e+03 7.376e+03
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 2.864e+01 1.939e+01 2.844e+01 3.869e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 1.116e+05 6.618e+04 1.082e+05 1.769e+05
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 5.730e+01 5.316e+01 5.730e+01 6.090e+01
#> parameters
#> s(age).tau21 5.464e+04
#> s(age).alpha NA
#> s(age).edf 5.881e+00
#> s(wk_ini).tau21 1.166e+05
#> s(wk_ini).alpha NA
#> s(wk_ini).edf 8.603e+00
#> s(rivwk_conception).tau21 3.230e-01
#> s(rivwk_conception).alpha NA
#> s(rivwk_conception).edf 2.063e+00
#> s(res_muni).tau21 3.729e+03
#> s(res_muni).alpha NA
#> s(res_muni).edf 4.240e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 1.090e+04
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha NA
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 3.926e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 8.470e+03
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha NA
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 3.971e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 3.186e+05
#> 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 6.377e+01
#> ---
#> Formula sigma:
#> ---
#> sigma ~ sex + born_race + birth_place + study_years + consult_num +
#> s(age) + s(wk_ini) + s(res_muni, bs = "re")
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) 6.350902 6.337247 6.350610 6.366726 0.168
#> sexfemale -0.044834 -0.049633 -0.044746 -0.040103 -0.045
#> born_raceindigenous -0.074073 -0.083292 -0.073897 -0.065064 -0.074
#> born_raceNA 0.003023 -0.029275 0.002383 0.035611 0.003
#> birth_placeanother health center 0.047802 0.003670 0.047814 0.090789 0.047
#> birth_placehome 0.015113 0.007216 0.015159 0.023077 0.015
#> birth_placeother 0.087372 0.058241 0.087892 0.115886 0.086
#> birth_placeNA 0.021559 -0.008705 0.021234 0.053281 0.022
#> study_years1 - 3 -0.005366 -0.019164 -0.005103 0.008075 -0.005
#> study_years4 - 7 -0.002014 -0.014726 -0.001820 0.010564 -0.002
#> study_years8 - 11 0.025644 0.012229 0.025984 0.038795 0.026
#> study_years> 12 0.045591 0.027834 0.045485 0.062166 0.046
#> study_yearsNA -0.007150 -0.039975 -0.006873 0.026541 -0.007
#> consult_num1 - 3 -0.075669 -0.086735 -0.075741 -0.063678 -0.076
#> consult_num4 - 6 -0.155835 -0.166889 -0.155809 -0.143868 -0.156
#> consult_num> 7 -0.212554 -0.224049 -0.212744 -0.200100 -0.212
#> consult_numNA -0.011255 -0.039495 -0.010867 0.016002 -0.012
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.9393 0.7144 0.9866 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5% parameters
#> s(age).tau21 0.046685 0.007912 0.033372 0.163191 0.031
#> s(age).alpha 0.975315 0.848687 0.996519 1.000000 NA
#> s(age).edf 4.995966 3.639184 4.980607 6.336074 4.920
#> s(wk_ini).tau21 0.043720 0.009475 0.033226 0.151527 0.075
#> s(wk_ini).alpha 0.980771 0.887841 0.996616 1.000000 NA
#> s(wk_ini).edf 7.199591 5.649132 7.268384 8.438409 7.946
#> s(res_muni).tau21 38.304091 23.828566 37.063513 59.660053 38.146
#> s(res_muni).alpha 0.939796 0.708036 0.992603 1.000000 NA
#> s(res_muni).edf 42.999866 42.999806 42.999869 42.999920 43.000
#> ---
#> Sampler summary:
#> -
#> runtime = 61126.91
#> ---
#> Optimizer summary:
#> -
#> AICc = 4433567 edf = 298.5487 logLik = -2216485
#> logPost = -2218409 nobs = 291479 runtime = 1284.172Parametric effects
par(mar = c(4, 4, 0.5, 0), mfrow = c(1, 2), cex.axis = 0.7)
plot2d_bamlss(model, bwdata_model, model = "mu", term = "remoteness", grid = 50,
FUN = c95)
plot2d_bamlss(model, bwdata_model, model = "mu", 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
#> 10.591 0.315 11.020