Growth model with t-distribution and no rainfall: 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-25-growth-no-rain-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 = 8000, 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")
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) 3102.283 3041.208 3103.395 3157.566 3134.267
#> sexfemale -101.012 -104.264 -100.962 -97.606 -101.046
#> born_raceindigenous -50.208 -55.586 -50.196 -44.773 -50.152
#> born_raceNA 4.859 -14.839 4.964 23.863 4.666
#> birth_placeanother health center 11.923 -15.863 11.941 40.659 11.663
#> birth_placehome -62.748 -68.342 -62.687 -57.363 -62.693
#> birth_placeother -71.160 -90.462 -71.209 -52.012 -71.156
#> birth_placeNA -102.697 -123.639 -102.379 -81.951 -102.492
#> gestation_weeks32 - 36 -246.873 -253.327 -246.916 -240.382 -247.049
#> gestation_weeks28 - 31 -853.165 -879.892 -853.390 -827.124 -853.838
#> gestation_weeks22 - 27 -1772.918 -1818.968 -1772.764 -1723.943 -1772.845
#> marital_statusmarried 30.936 25.668 30.953 36.364 30.989
#> marital_statuswidow -11.533 -63.998 -11.915 39.861 -12.393
#> marital_statusdivorced -32.859 -79.202 -31.474 12.632 -32.821
#> marital_statusconsensual union 24.213 19.819 24.230 28.767 24.227
#> marital_statusNA 22.372 8.155 22.190 38.418 22.043
#> study_years1 - 3 60.665 51.959 60.824 68.938 60.675
#> study_years4 - 7 81.690 74.028 81.658 89.794 81.721
#> study_years8 - 11 75.396 67.095 75.211 84.003 75.487
#> study_years> 12 72.494 62.341 72.384 83.240 72.660
#> study_yearsNA 47.010 25.328 47.127 68.787 46.973
#> consult_num1 - 3 66.813 58.554 66.920 74.139 66.866
#> consult_num4 - 6 101.831 93.998 101.959 108.926 101.974
#> consult_num> 7 146.946 138.967 147.016 154.609 147.100
#> consult_numNA 21.169 2.092 21.141 41.225 21.102
#> remoteness -4.895 -66.232 -5.968 53.484 -10.860
#> prop_tap_toilet 132.840 17.150 130.836 264.759 -70.476
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8290 0.2997 0.9297 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5% parameters
#> s(age).tau21 1.644e+05 4.203e+04 1.282e+05 5.043e+05 4.309e+04
#> s(age).alpha 9.809e-01 8.765e-01 9.995e-01 1.000e+00 NA
#> s(age).edf 6.516e+00 5.780e+00 6.519e+00 7.329e+00 5.784e+00
#> s(wk_ini).tau21 9.611e+04 2.649e+04 7.651e+04 2.897e+05 2.016e+05
#> s(wk_ini).alpha 9.966e-01 9.810e-01 9.999e-01 1.000e+00 NA
#> s(wk_ini).edf 8.461e+00 7.830e+00 8.501e+00 8.867e+00 8.787e+00
#> s(rivwk_conception).tau21 2.232e+00 1.887e-01 1.363e+00 1.039e+01 1.909e+00
#> s(rivwk_conception).alpha 9.991e-01 9.946e-01 1.000e+00 1.000e+00 NA
#> s(rivwk_conception).edf 3.481e+00 1.808e+00 3.437e+00 5.608e+00 3.749e+00
#> s(res_muni).tau21 1.908e+03 1.228e+03 1.841e+03 2.884e+03 4.152e+03
#> s(res_muni).alpha 9.857e-01 9.307e-01 9.999e-01 1.000e+00 NA
#> s(res_muni).edf 4.192e+01 4.142e+01 4.193e+01 4.232e+01 4.252e+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.1077914 6.0893580 6.1077604 6.1272758 0.082
#> sexfemale -0.0410126 -0.0463576 -0.0410585 -0.0352277 -0.041
#> born_raceindigenous -0.0566448 -0.0664529 -0.0566675 -0.0461373 -0.057
#> born_raceNA -0.0008566 -0.0385714 -0.0008478 0.0346444 -0.001
#> birth_placeanother health center 0.0630262 0.0114596 0.0632751 0.1117714 0.063
#> birth_placehome 0.0522655 0.0421192 0.0523469 0.0623431 0.052
#> birth_placeother 0.0610985 0.0267787 0.0612330 0.0954881 0.061
#> birth_placeNA 0.0869311 0.0513823 0.0865297 0.1217677 0.087
#> study_years1 - 3 0.0026571 -0.0131167 0.0026073 0.0183816 0.002
#> study_years4 - 7 0.0044598 -0.0111246 0.0042728 0.0188575 0.004
#> study_years8 - 11 0.0149874 -0.0011356 0.0150063 0.0308461 0.015
#> study_years> 12 0.0239234 0.0034686 0.0237533 0.0446986 0.024
#> study_yearsNA -0.0222662 -0.0611596 -0.0229173 0.0162926 -0.023
#> consult_num1 - 3 -0.0385886 -0.0520717 -0.0387949 -0.0246325 -0.038
#> consult_num4 - 6 -0.0803982 -0.0930151 -0.0803824 -0.0670049 -0.080
#> consult_num> 7 -0.1071433 -0.1215766 -0.1071288 -0.0933517 -0.107
#> consult_numNA 0.0032929 -0.0318429 0.0034741 0.0375507 0.002
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.9628 0.8114 0.9932 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5% parameters
#> s(age).tau21 9.550e-03 2.378e-04 4.198e-03 5.701e-02 0.005
#> s(age).alpha 9.922e-01 9.406e-01 9.997e-01 1.000e+00 NA
#> s(age).edf 3.054e+00 1.409e+00 2.994e+00 5.179e+00 3.141
#> s(wk_ini).tau21 1.318e-02 7.317e-04 8.639e-03 5.202e-02 0.018
#> s(wk_ini).alpha 9.959e-01 9.749e-01 9.999e-01 1.000e+00 NA
#> s(wk_ini).edf 5.290e+00 2.991e+00 5.344e+00 7.338e+00 6.173
#> s(res_muni).tau21 3.467e+01 2.243e+01 3.403e+01 5.190e+01 36.211
#> s(res_muni).alpha 9.359e-01 6.752e-01 9.982e-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.117 2.087 2.117 2.144 2.121
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8851 0.3133 0.9895 1
#> ---
#> Sampler summary:
#> -
#> runtime = 52206.79
#> ---
#> Optimizer summary:
#> -
#> AICc = 4408851 edf = 158.1527 logLik = -2204267
#> logPost = -2205082 nobs = 291479 runtime = 286.17Parametric 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.762 0.319 10.166