Growth model without 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, "lbw-25-growth-no-rain.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 = "binomial", data = bwdata_model,
#> cores = 4, combine = FALSE, light = TRUE, n.iter = 8000,
#> burnin = 0)
#> ---
#> Family: binomial
#> Link function: pi = logit
#> *---
#> Formula pi:
#> ---
#> lbw ~ 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) -2.340651 -2.634095 -2.340099 -2.036580 0.522
#> sexfemale 0.237686 0.206610 0.237310 0.268643 0.238
#> born_raceindigenous 0.019058 -0.031709 0.019045 0.071749 0.021
#> born_raceNA -0.017311 -0.207995 -0.019638 0.188155 -0.013
#> birth_placeanother health center -0.027217 -0.285194 -0.029861 0.234747 -0.022
#> birth_placehome -0.008129 -0.056217 -0.008223 0.038796 -0.006
#> birth_placeother 0.203643 0.053623 0.203288 0.359849 0.205
#> birth_placeNA 0.024601 -0.146251 0.025108 0.193526 0.023
#> gestation_weeks32 - 36 1.560069 1.520436 1.560021 1.599138 1.563
#> gestation_weeks28 - 31 3.120446 3.022957 3.119627 3.217105 3.125
#> gestation_weeks22 - 27 3.858972 3.684736 3.859895 4.025106 3.862
#> marital_statusmarried -0.162582 -0.215580 -0.162589 -0.106140 -0.164
#> marital_statuswidow -0.316044 -0.860433 -0.321352 0.260190 -0.325
#> marital_statusdivorced -0.103366 -0.569419 -0.101525 0.354449 -0.112
#> marital_statusconsensual union -0.140115 -0.184294 -0.139768 -0.096036 -0.141
#> marital_statusNA 0.020092 -0.130347 0.018642 0.161176 0.028
#> study_years1 - 3 -0.247338 -0.315837 -0.246816 -0.179326 -0.247
#> study_years4 - 7 -0.300448 -0.368752 -0.299326 -0.233219 -0.301
#> study_years8 - 11 -0.248919 -0.317755 -0.249286 -0.180925 -0.252
#> study_years> 12 -0.087228 -0.183909 -0.086602 0.006483 -0.092
#> study_yearsNA -0.188997 -0.388407 -0.189301 0.001532 -0.180
#> consult_num1 - 3 -0.436384 -0.500191 -0.435827 -0.374882 -0.435
#> consult_num4 - 6 -0.642464 -0.707034 -0.641543 -0.584040 -0.641
#> consult_num> 7 -0.912524 -0.980567 -0.911757 -0.846120 -0.911
#> consult_numNA -0.045159 -0.190640 -0.044193 0.099581 -0.044
#> remoteness 0.173398 -0.123179 0.176617 0.491875 -0.122
#> prop_tap_toilet -0.280567 -0.880353 -0.279022 0.280466 -0.142
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.7565 0.0843 0.8599 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5% parameters
#> s(age).tau21 3.398e+00 8.461e-01 2.505e+00 1.173e+01 0.666
#> s(age).alpha 8.935e-01 3.088e-01 9.858e-01 1.000e+00 NA
#> s(age).edf 5.693e+00 4.728e+00 5.673e+00 6.810e+00 4.521
#> s(wk_ini).tau21 3.988e+00 8.367e-01 3.122e+00 1.218e+01 9.025
#> s(wk_ini).alpha 9.718e-01 8.399e-01 9.989e-01 1.000e+00 NA
#> s(wk_ini).edf 7.954e+00 6.782e+00 8.033e+00 8.690e+00 8.585
#> s(rivwk_conception).tau21 2.098e-04 3.489e-05 1.361e-04 8.544e-04 0.000
#> s(rivwk_conception).alpha 9.914e-01 9.427e-01 9.998e-01 1.000e+00 NA
#> s(rivwk_conception).edf 3.585e+00 2.275e+00 3.489e+00 5.453e+00 2.298
#> s(res_muni).tau21 3.521e-02 2.107e-02 3.405e-02 5.646e-02 8.080
#> s(res_muni).alpha 7.725e-01 2.395e-01 8.415e-01 1.000e+00 NA
#> s(res_muni).edf 3.799e+01 3.571e+01 3.804e+01 3.987e+01 42.975
#> ---
#> Sampler summary:
#> -
#> runtime = 10518.15
#> ---
#> Optimizer summary:
#> -
#> AICc = 120471.2 edf = 85.3796 logLik = -60150.2
#> logPost = -60436.25 nobs = 291479 runtime = 35.215Parametric effects
par(mar = c(4, 4, 0.5, 0), mfrow = c(1, 2), cex.axis = 0.7)
plot2d_bamlss(model, bwdata_model, model = "pi", term = "remoteness", grid = 50,
FUN = c95)
plot2d_bamlss(model, bwdata_model, model = "pi", 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
#> 6.104 0.211 6.333