Growth 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-20-growth-re-t.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 = TF, data = bwdata_model, cores = 4,
#> combine = FALSE, light = TRUE, n.iter = 30, burnin = 0)
#> ---
#> Family: TF
#> Link function: mu = identity, sigma = log, nu = log
#> *---
#> Formula mu:
#> ---
#> born_weight ~ sex + born_race + birth_place + gestation_weeks +
#> 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) 3104.0220 3045.7682 3105.0644 3157.4153 3133.366
#> sexfemale -101.1032 -104.4334 -101.0706 -97.8916 -101.071
#> born_raceindigenous -49.7906 -55.2326 -49.8641 -43.9264 -49.465
#> born_raceNA 3.2194 -17.3920 3.3628 23.6349 2.563
#> birth_placeanother health center 12.3545 -17.9212 12.5570 41.9300 12.979
#> birth_placehome -63.1751 -68.4246 -63.1610 -57.6719 -63.177
#> birth_placeother -72.2150 -91.0084 -72.5160 -52.4770 -72.581
#> birth_placeNA -103.9562 -125.0113 -103.2277 -83.6193 -102.641
#> gestation_weeks32 - 36 -247.3127 -253.8111 -247.2296 -240.9529 -247.528
#> gestation_weeks28 - 31 -853.2885 -880.5655 -852.9429 -826.9516 -853.824
#> gestation_weeks22 - 27 -1774.3791 -1819.6321 -1775.1359 -1728.3186 -1773.369
#> marital_statusmarried 31.4146 25.7453 31.4685 36.9698 31.559
#> marital_statuswidow -12.3297 -66.6040 -12.4580 41.5353 -12.517
#> marital_statusdivorced -30.8328 -74.6428 -31.1723 16.6780 -30.619
#> marital_statusconsensual union 24.9866 20.2013 24.9766 29.5447 25.270
#> marital_statusNA 24.2958 9.0298 24.1021 40.0293 24.485
#> study_years1 - 3 60.7484 51.7162 60.7632 69.4668 60.715
#> study_years4 - 7 81.5980 73.2929 81.7443 89.4107 81.587
#> study_years8 - 11 75.2166 66.8088 75.3807 83.7145 75.328
#> study_years> 12 72.9782 61.7735 72.8941 84.0968 73.293
#> study_yearsNA 47.0688 25.4032 47.0085 68.0343 46.322
#> consult_num1 - 3 66.5303 58.0885 66.5804 74.1679 66.913
#> consult_num4 - 6 101.3510 93.6493 101.2517 109.0798 101.757
#> consult_num> 7 146.7547 138.3480 146.8090 154.4174 147.280
#> consult_numNA 20.6684 0.2648 20.9783 41.3571 20.672
#> remoteness -4.5862 -64.4733 -3.5486 58.3201 -10.486
#> prop_tap_toilet 126.6656 9.5268 124.6175 244.1859 -69.015
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8349 0.2952 0.9368 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 1.699e+05 4.236e+04 1.287e+05 5.279e+05
#> s(age).alpha 9.770e-01 8.529e-01 9.977e-01 1.000e+00
#> s(age).edf 6.551e+00 5.803e+00 6.556e+00 7.363e+00
#> s(wk_ini).tau21 1.144e+05 2.700e+04 9.055e+04 3.383e+05
#> s(wk_ini).alpha 9.967e-01 9.823e-01 9.999e-01 1.000e+00
#> s(wk_ini).edf 8.511e+00 7.879e+00 8.556e+00 8.870e+00
#> s(rivwk_conception).tau21 2.297e+00 2.246e-01 1.386e+00 9.826e+00
#> s(rivwk_conception).alpha 9.991e-01 9.950e-01 1.000e+00 1.000e+00
#> s(rivwk_conception).edf 3.563e+00 1.948e+00 3.480e+00 5.749e+00
#> s(res_muni).tau21 1.900e+03 1.210e+03 1.836e+03 2.910e+03
#> s(res_muni).alpha 9.852e-01 9.255e-01 1.000e+00 1.000e+00
#> s(res_muni).edf 4.191e+01 4.141e+01 4.193e+01 4.231e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 1.053e+03 8.708e-01 7.538e+02 3.669e+03
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 9.920e-01 9.505e-01 1.000e+00 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 1.518e+01 2.061e+00 1.523e+01 2.806e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 9.609e+02 2.606e+02 8.596e+02 2.347e+03
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 9.916e-01 9.546e-01 9.996e-01 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 1.848e+01 1.117e+01 1.808e+01 2.693e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 5.352e+03 9.113e+02 4.880e+03 1.309e+04
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 9.753e-01 8.760e-01 9.977e-01 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 2.809e+01 1.661e+01 2.853e+01 3.793e+01
#> parameters
#> s(age).tau21 4.398e+04
#> s(age).alpha NA
#> s(age).edf 5.800e+00
#> s(wk_ini).tau21 2.252e+05
#> s(wk_ini).alpha NA
#> s(wk_ini).edf 8.808e+00
#> s(rivwk_conception).tau21 1.881e+00
#> s(rivwk_conception).alpha NA
#> s(rivwk_conception).edf 3.735e+00
#> s(res_muni).tau21 4.034e+03
#> s(res_muni).alpha NA
#> s(res_muni).edf 4.251e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 4.938e+03
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha NA
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 3.164e+01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 9.154e+02
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha NA
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 1.871e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 2.566e+04
#> 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 4.504e+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.1083187 6.0884242 6.1085467 6.1263361 0.083
#> sexfemale -0.0409832 -0.0469649 -0.0409314 -0.0352805 -0.041
#> born_raceindigenous -0.0571824 -0.0672195 -0.0574214 -0.0466511 -0.057
#> born_raceNA -0.0006248 -0.0365270 -0.0008937 0.0357110 -0.001
#> birth_placeanother health center 0.0630651 0.0133303 0.0632764 0.1155492 0.063
#> birth_placehome 0.0522271 0.0425654 0.0521730 0.0618733 0.052
#> birth_placeother 0.0611747 0.0264977 0.0609290 0.0956036 0.061
#> birth_placeNA 0.0866361 0.0516586 0.0862344 0.1216696 0.087
#> study_years1 - 3 0.0025600 -0.0129378 0.0026473 0.0185059 0.002
#> study_years4 - 7 0.0039396 -0.0104395 0.0039979 0.0180910 0.004
#> study_years8 - 11 0.0145051 -0.0007927 0.0147907 0.0291505 0.015
#> study_years> 12 0.0230982 0.0041006 0.0232580 0.0419995 0.024
#> study_yearsNA -0.0231307 -0.0628983 -0.0227722 0.0148448 -0.022
#> consult_num1 - 3 -0.0390029 -0.0535609 -0.0387638 -0.0257134 -0.039
#> consult_num4 - 6 -0.0805355 -0.0936372 -0.0804452 -0.0672630 -0.081
#> consult_num> 7 -0.1073855 -0.1213685 -0.1071252 -0.0940621 -0.107
#> consult_numNA 0.0012127 -0.0334769 0.0016119 0.0356728 0.001
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.9612 0.7845 0.9997 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5% parameters
#> s(age).tau21 1.025e-02 2.406e-04 3.927e-03 6.026e-02 0.005
#> s(age).alpha 9.920e-01 9.360e-01 9.998e-01 1.000e+00 NA
#> s(age).edf 3.040e+00 1.412e+00 2.938e+00 5.226e+00 3.134
#> s(wk_ini).tau21 1.373e-02 8.919e-04 8.595e-03 5.393e-02 0.018
#> s(wk_ini).alpha 9.957e-01 9.731e-01 1.000e+00 1.000e+00 NA
#> s(wk_ini).edf 5.362e+00 3.215e+00 5.413e+00 7.308e+00 6.192
#> s(res_muni).tau21 3.758e+01 2.401e+01 3.671e+01 5.609e+01 36.266
#> s(res_muni).alpha 9.368e-01 7.131e-01 9.932e-01 1.000e+00 NA
#> s(res_muni).edf 4.300e+01 4.300e+01 4.300e+01 4.300e+01 43.000
#> ---
#> Formula nu:
#> ---
#> nu ~ 1
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) 2.116 2.090 2.115 2.146 2.121
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8788 0.3704 0.9826 1
#> ---
#> Sampler summary:
#> -
#> runtime = 369.158
#> ---
#> Optimizer summary:
#> -
#> AICc = 4408738 edf = 253.5741 logLik = -2204115
#> logPost = -2205868 nobs = 291479 runtime = 846.639Parametric 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.532 0.342 10.944