Full 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-15-full-no-rain-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 = 8000, 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")
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) 3095.420 3033.270 3095.102 3159.944 3127.568
#> sexfemale -100.281 -103.773 -100.274 -97.030 -100.216
#> born_raceindigenous -58.568 -64.121 -58.443 -53.162 -58.511
#> born_raceNA 4.681 -15.582 4.545 26.025 4.359
#> birth_placeanother health center 5.488 -25.888 5.378 34.775 5.119
#> birth_placehome -73.592 -79.091 -73.582 -68.363 -73.643
#> birth_placeother -86.975 -106.597 -87.207 -68.001 -87.222
#> birth_placeNA -114.220 -133.817 -114.288 -93.595 -114.851
#> marital_statusmarried 30.833 25.128 30.792 36.539 30.892
#> marital_statuswidow -11.608 -62.532 -12.545 41.436 -11.513
#> marital_statusdivorced -29.539 -73.188 -30.041 19.650 -29.778
#> marital_statusconsensual union 25.598 21.146 25.619 30.144 25.322
#> marital_statusNA 22.139 6.658 21.976 38.575 21.761
#> study_years1 - 3 62.530 54.099 62.686 71.076 62.663
#> study_years4 - 7 84.173 75.517 84.303 91.602 84.192
#> study_years8 - 11 77.116 68.243 77.227 85.325 77.226
#> study_years> 12 75.061 64.154 75.033 85.708 75.330
#> study_yearsNA 52.047 31.220 52.066 73.571 52.216
#> consult_num1 - 3 49.668 41.287 49.677 57.982 49.707
#> consult_num4 - 6 92.068 84.414 92.082 99.774 92.112
#> consult_num> 7 150.903 142.301 150.857 159.471 150.801
#> consult_numNA 17.450 -3.412 17.611 37.651 17.398
#> remoteness -7.630 -77.072 -9.474 68.796 -18.465
#> prop_tap_toilet 111.629 -13.420 110.532 239.062 -73.364
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.9736 0.8703 0.9969 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5% parameters
#> s(age).tau21 1.720e+05 4.808e+04 1.330e+05 5.003e+05 3.905e+04
#> s(age).alpha 9.752e-01 8.366e-01 9.980e-01 1.000e+00 NA
#> s(age).edf 6.572e+00 5.823e+00 6.553e+00 7.403e+00 5.719e+00
#> s(wk_ini).tau21 2.170e+04 6.258e+02 1.413e+04 9.012e+04 1.119e+05
#> s(wk_ini).alpha 9.970e-01 9.841e-01 9.997e-01 1.000e+00 NA
#> s(wk_ini).edf 6.932e+00 3.933e+00 7.264e+00 8.559e+00 8.630e+00
#> s(rivwk_conception).tau21 7.778e-01 3.876e-03 4.156e-01 3.695e+00 3.130e-01
#> s(rivwk_conception).alpha 9.993e-01 9.952e-01 1.000e+00 1.000e+00 NA
#> s(rivwk_conception).edf 2.390e+00 1.275e-01 2.361e+00 4.476e+00 2.131e+00
#> s(res_muni).tau21 1.909e+03 1.218e+03 1.860e+03 2.994e+03 3.582e+03
#> s(res_muni).alpha 9.820e-01 9.116e-01 1.000e+00 1.000e+00 NA
#> s(res_muni).edf 4.190e+01 4.141e+01 4.191e+01 4.234e+01 4.244e+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.1211668 6.1013504 6.1212178 6.1410708 0.113
#> sexfemale -0.0433806 -0.0497822 -0.0432779 -0.0375054 -0.043
#> born_raceindigenous -0.0704755 -0.0810196 -0.0703877 -0.0594267 -0.070
#> born_raceNA -0.0069248 -0.0454126 -0.0074375 0.0300311 -0.006
#> birth_placeanother health center 0.0438438 -0.0119166 0.0439106 0.0949864 0.043
#> birth_placehome 0.0253905 0.0152728 0.0253180 0.0358169 0.025
#> birth_placeother 0.0704168 0.0361950 0.0704916 0.1060995 0.069
#> birth_placeNA 0.0242636 -0.0140943 0.0238164 0.0652116 0.025
#> study_years1 - 3 0.0003537 -0.0162798 0.0003427 0.0166258 0.001
#> study_years4 - 7 0.0043230 -0.0095131 0.0040975 0.0192806 0.005
#> study_years8 - 11 0.0224433 0.0063708 0.0223101 0.0387103 0.023
#> study_years> 12 0.0389814 0.0186381 0.0387643 0.0589114 0.039
#> study_yearsNA -0.0094619 -0.0512963 -0.0093080 0.0290985 -0.010
#> consult_num1 - 3 -0.0537736 -0.0685076 -0.0538481 -0.0377746 -0.054
#> consult_num4 - 6 -0.1063576 -0.1200019 -0.1064995 -0.0913345 -0.107
#> consult_num> 7 -0.1434501 -0.1576437 -0.1436210 -0.1282759 -0.144
#> consult_numNA -0.0135890 -0.0490329 -0.0128005 0.0214813 -0.014
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.9577 0.7821 0.9990 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5% parameters
#> s(age).tau21 0.030806 0.003449 0.018107 0.153861 0.022
#> s(age).alpha 0.987259 0.913283 0.998883 1.000000 NA
#> s(age).edf 4.162761 2.760010 4.146511 5.802593 4.264
#> s(wk_ini).tau21 0.021550 0.001923 0.013497 0.093290 0.035
#> s(wk_ini).alpha 0.993094 0.952078 0.999841 1.000000 NA
#> s(wk_ini).edf 5.796022 3.731669 5.853465 7.746031 6.879
#> s(res_muni).tau21 33.753080 22.164898 32.608406 51.161699 36.364
#> s(res_muni).alpha 0.931994 0.700884 0.988449 1.000000 NA
#> s(res_muni).edf 42.999777 42.999674 42.999781 42.999861 43.000
#> ---
#> Formula nu:
#> ---
#> nu ~ 1
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) 1.876 1.854 1.876 1.898 1.88
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8848 0.3869 0.9877 1
#> ---
#> Sampler summary:
#> -
#> runtime = 52095.72
#> ---
#> Optimizer summary:
#> -
#> AICc = 4420540 edf = 155.0627 logLik = -2210115
#> logPost = -2210901 nobs = 291479 runtime = 589.692Parametric 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.390 0.297 9.774