Full 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, "pre-15-full-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:
#> ---
#> premature ~ 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) -3.039e+00 -3.553e+00 -3.052e+00 -2.457e+00 -0.168
#> sexfemale -3.524e-02 -6.012e-02 -3.539e-02 -1.013e-02 -0.035
#> born_raceindigenous 2.606e-01 2.209e-01 2.605e-01 3.022e-01 0.259
#> born_raceNA 5.436e-02 -7.917e-02 5.115e-02 2.001e-01 0.053
#> birth_placeanother health center 1.594e-01 -2.275e-02 1.566e-01 3.661e-01 0.169
#> birth_placehome 2.048e-01 1.653e-01 2.045e-01 2.448e-01 0.204
#> birth_placeother 3.151e-01 1.922e-01 3.174e-01 4.245e-01 0.314
#> birth_placeNA -3.727e-05 -1.228e-01 7.218e-05 1.201e-01 0.004
#> marital_statusmarried -3.043e-02 -8.118e-02 -3.026e-02 1.461e-02 -0.031
#> marital_statuswidow 3.194e-02 -3.939e-01 4.265e-02 4.820e-01 0.036
#> marital_statusdivorced -2.153e-01 -5.956e-01 -2.245e-01 2.191e-01 -0.226
#> marital_statusconsensual union -5.406e-02 -8.689e-02 -5.296e-02 -2.440e-02 -0.056
#> marital_statusNA -1.059e-02 -1.321e-01 -1.027e-02 1.126e-01 -0.014
#> study_years1 - 3 -6.971e-02 -1.371e-01 -6.992e-02 -2.596e-03 -0.071
#> study_years4 - 7 -7.154e-02 -1.298e-01 -7.065e-02 -1.281e-02 -0.073
#> study_years8 - 11 9.068e-03 -5.346e-02 9.185e-03 7.124e-02 0.005
#> study_years> 12 4.044e-02 -5.432e-02 4.088e-02 1.281e-01 0.036
#> study_yearsNA -1.806e-01 -3.202e-01 -1.779e-01 -4.086e-02 -0.188
#> consult_num1 - 3 5.023e-01 4.390e-01 5.027e-01 5.633e-01 0.501
#> consult_num4 - 6 2.254e-01 1.659e-01 2.247e-01 2.875e-01 0.224
#> consult_num> 7 -5.600e-01 -6.273e-01 -5.605e-01 -4.932e-01 -0.562
#> consult_numNA 1.207e-01 -2.526e-02 1.175e-01 2.809e-01 0.119
#> remoteness 2.903e-02 -5.176e-01 3.370e-02 5.782e-01 -0.152
#> prop_tap_toilet 9.300e-01 -2.311e-01 9.576e-01 1.908e+00 0.251
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.8169 0.2048 0.9084 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5% parameters
#> s(age).tau21 2.051e+00 4.559e-01 1.488e+00 6.782e+00 0.447
#> s(age).alpha 9.301e-01 4.974e-01 9.912e-01 1.000e+00 NA
#> s(age).edf 5.552e+00 4.495e+00 5.548e+00 6.661e+00 4.462
#> s(wk_ini).tau21 3.623e+01 1.232e+01 2.874e+01 1.131e+02 27.363
#> s(wk_ini).alpha 9.692e-01 8.241e-01 9.929e-01 1.000e+00 NA
#> s(wk_ini).edf 8.869e+00 8.724e+00 8.879e+00 8.962e+00 8.867
#> s(rivwk_conception).tau21 1.169e-03 3.129e-04 9.011e-04 3.684e-03 0.001
#> s(rivwk_conception).alpha 9.864e-01 9.188e-01 9.989e-01 1.000e+00 NA
#> s(rivwk_conception).edf 5.992e+00 4.802e+00 5.988e+00 7.156e+00 5.387
#> s(res_muni).tau21 1.192e-01 7.366e-02 1.151e-01 1.872e-01 7.258
#> s(res_muni).alpha 7.810e-01 2.323e-01 8.764e-01 1.000e+00 NA
#> s(res_muni).edf 4.180e+01 4.121e+01 4.182e+01 4.226e+01 42.980
#> ---
#> Sampler summary:
#> -
#> runtime = 20560.72
#> ---
#> Optimizer summary:
#> -
#> AICc = 160987.1 edf = 85.6955 logLik = -80407.83
#> logPost = -80685.28 nobs = 291479 runtime = 170.599Parametric 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, contour = TRUE, image = TRUE)
Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 6.093 0.203 6.312