Full model with t-distribution: 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-10-full-re-t.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, family = TF, data = bwdata_model, cores = 4,
#> combine = FALSE, light = TRUE, n.iter = 1000, burnin = 0)
#> ---
#> Family: TF
#> Link function: mu = identity, sigma = log, nu = 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) 3098.171 3040.026 3097.444 3157.314 3127.974
#> sexfemale -100.385 -103.728 -100.400 -97.037 -100.296
#> born_raceindigenous -57.535 -62.931 -57.499 -52.200 -57.616
#> born_raceNA 4.253 -16.417 4.248 25.132 3.908
#> birth_placeanother health center 6.086 -24.108 6.249 36.214 6.186
#> birth_placehome -74.150 -79.353 -74.109 -68.995 -73.976
#> birth_placeother -87.903 -108.560 -87.588 -66.551 -88.372
#> birth_placeNA -114.229 -133.157 -114.416 -95.789 -113.776
#> marital_statusmarried 31.571 26.119 31.549 37.191 31.602
#> marital_statuswidow -10.857 -59.360 -12.091 41.909 -10.646
#> marital_statusdivorced -28.199 -73.846 -28.416 16.988 -28.633
#> marital_statusconsensual union 26.635 22.169 26.712 31.098 26.718
#> marital_statusNA 23.514 7.630 23.654 39.569 23.824
#> study_years1 - 3 63.025 54.600 63.152 71.531 63.155
#> study_years4 - 7 84.305 76.040 84.312 92.507 84.357
#> study_years8 - 11 77.240 68.753 77.242 85.401 77.242
#> study_years> 12 75.968 65.154 75.837 87.558 76.017
#> study_yearsNA 51.053 29.224 51.322 71.362 50.534
#> consult_num1 - 3 49.722 42.390 49.681 57.625 49.774
#> consult_num4 - 6 91.801 84.511 91.902 98.832 91.803
#> consult_num> 7 150.700 142.876 150.730 158.690 150.667
#> consult_numNA 16.570 -4.741 16.667 36.959 16.661
#> remoteness -12.001 -72.632 -12.886 49.346 -19.643
#> prop_tap_toilet 106.267 -16.694 106.749 226.257 -75.855
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.9735 0.8658 0.9979 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 1.812e+05 5.026e+04 1.389e+05 6.040e+05
#> s(age).alpha 9.762e-01 8.452e-01 9.992e-01 1.000e+00
#> s(age).edf 6.619e+00 5.917e+00 6.604e+00 7.400e+00
#> s(wk_ini).tau21 7.258e+03 1.651e+02 2.413e+03 4.212e+04
#> s(wk_ini).alpha 9.979e-01 9.873e-01 1.000e+00 1.000e+00
#> s(wk_ini).edf 5.400e+00 2.885e+00 5.207e+00 8.241e+00
#> s(rivwk_conception).tau21 8.820e-01 4.163e-03 4.466e-01 4.319e+00
#> s(rivwk_conception).alpha 9.993e-01 9.950e-01 1.000e+00 1.000e+00
#> s(rivwk_conception).edf 2.445e+00 1.767e-01 2.414e+00 4.572e+00
#> s(res_muni).tau21 1.823e+03 1.172e+03 1.748e+03 2.845e+03
#> s(res_muni).alpha 9.822e-01 9.111e-01 9.978e-01 1.000e+00
#> s(res_muni).edf 4.186e+01 4.129e+01 4.188e+01 4.230e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 5.197e+02 7.662e-01 3.501e+02 2.159e+03
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 9.926e-01 9.594e-01 9.990e-01 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 1.083e+01 2.087e+00 1.063e+01 2.282e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 1.543e+03 4.377e+02 1.355e+03 3.559e+03
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 9.880e-01 9.382e-01 9.998e-01 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 2.207e+01 1.398e+01 2.191e+01 3.117e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 5.760e+04 2.945e+04 5.605e+04 9.661e+04
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 9.256e-01 6.866e-01 9.874e-01 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 5.229e+01 4.645e+01 5.244e+01 5.739e+01
#> parameters
#> s(age).tau21 3.874e+04
#> s(age).alpha NA
#> s(age).edf 5.714e+00
#> s(wk_ini).tau21 1.079e+05
#> s(wk_ini).alpha NA
#> s(wk_ini).edf 8.619e+00
#> s(rivwk_conception).tau21 3.250e-01
#> s(rivwk_conception).alpha NA
#> s(rivwk_conception).edf 2.160e+00
#> s(res_muni).tau21 3.514e+03
#> s(res_muni).alpha NA
#> s(res_muni).edf 4.243e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 7.230e+02
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha NA
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 1.469e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 2.620e+03
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha NA
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 2.800e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 1.834e+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.163e+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.1216358 6.1023620 6.1217264 6.1410523 0.112
#> sexfemale -0.0433255 -0.0492370 -0.0433804 -0.0372285 -0.043
#> born_raceindigenous -0.0705560 -0.0814157 -0.0706281 -0.0597724 -0.071
#> born_raceNA -0.0064638 -0.0446406 -0.0060706 0.0340295 -0.006
#> birth_placeanother health center 0.0434365 -0.0047400 0.0443900 0.0977410 0.043
#> birth_placehome 0.0255851 0.0154573 0.0255527 0.0352919 0.026
#> birth_placeother 0.0715635 0.0380244 0.0710848 0.1079411 0.071
#> birth_placeNA 0.0244067 -0.0098765 0.0245579 0.0618413 0.025
#> study_years1 - 3 0.0002604 -0.0163886 0.0005253 0.0157620 0.000
#> study_years4 - 7 0.0039563 -0.0115420 0.0041544 0.0186836 0.004
#> study_years8 - 11 0.0218914 0.0053769 0.0218823 0.0374032 0.022
#> study_years> 12 0.0386046 0.0193703 0.0382010 0.0592433 0.039
#> study_yearsNA -0.0099752 -0.0474218 -0.0103115 0.0281008 -0.010
#> consult_num1 - 3 -0.0544704 -0.0686053 -0.0546685 -0.0398027 -0.054
#> consult_num4 - 6 -0.1068352 -0.1205922 -0.1067169 -0.0920768 -0.107
#> consult_num> 7 -0.1436410 -0.1582454 -0.1437191 -0.1278226 -0.144
#> consult_numNA -0.0151018 -0.0504646 -0.0155573 0.0189507 -0.015
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.9571 0.7741 0.9985 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5% parameters
#> s(age).tau21 0.035074 0.003317 0.020955 0.140782 0.023
#> s(age).alpha 0.987775 0.903588 1.000000 1.000000 NA
#> s(age).edf 4.251902 2.802336 4.201711 5.931355 4.314
#> s(wk_ini).tau21 0.021786 0.002221 0.015113 0.079930 0.035
#> s(wk_ini).alpha 0.993307 0.960458 1.000000 1.000000 NA
#> s(wk_ini).edf 5.894105 3.965739 5.921776 7.738381 6.859
#> s(res_muni).tau21 33.630737 21.538703 32.615577 52.416922 36.391
#> s(res_muni).alpha 0.930053 0.678384 0.988716 1.000000 NA
#> s(res_muni).edf 42.999774 42.999655 42.999779 42.999860 43.000
#> ---
#> Formula nu:
#> ---
#> nu ~ 1
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) 1.877 1.854 1.877 1.899 1.881
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8838 0.3908 0.9789 1
#> ---
#> Sampler summary:
#> -
#> runtime = 27704.27
#> ---
#> Optimizer summary:
#> -
#> AICc = 4420293 edf = 259.4192 logLik = -2209887
#> logPost = -2211647 nobs = 291479 runtime = 5561.42Parametric 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
#> 9.987 0.338 10.359