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.construct
library(gamlss.dist)
#> Loading required package: MASS
path_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.215

Parametric 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