Growth model: 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-20-growth-re.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 = 7000,
#> 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") +
#> s(neg_mbsi_mean_1wk, pos_mbsi_mean_1wk, k = 70) + s(neg_mckee_mean_8wk,
#> pos_mckee_mean_8wk, k = 70) + s(neg_ext_mbsi_mean_8wk, pos_ext_mbsi_mean_8wk,
#> k = 70)
#> -
#> Parametric coefficients:
#> Mean 2.5% 50% 97.5% parameters
#> (Intercept) -2.351877 -2.638702 -2.353068 -2.074827 0.480
#> sexfemale 0.238482 0.206895 0.238678 0.271451 0.239
#> born_raceindigenous 0.017756 -0.032343 0.017957 0.069437 0.021
#> born_raceNA -0.014844 -0.212361 -0.016477 0.184620 -0.008
#> birth_placeanother health center -0.021987 -0.262325 -0.020070 0.225628 -0.022
#> birth_placehome -0.005309 -0.051500 -0.005538 0.042714 -0.002
#> birth_placeother 0.208708 0.043737 0.208187 0.367929 0.213
#> birth_placeNA 0.020601 -0.143766 0.020962 0.178328 0.023
#> gestation_weeks32 - 36 1.564372 1.523728 1.564473 1.603788 1.567
#> gestation_weeks28 - 31 3.127436 3.032879 3.128571 3.223497 3.130
#> gestation_weeks22 - 27 3.858503 3.686736 3.857013 4.024766 3.866
#> marital_statusmarried -0.160349 -0.213905 -0.160619 -0.104303 -0.164
#> marital_statuswidow -0.300300 -0.861132 -0.318423 0.307211 -0.321
#> marital_statusdivorced -0.098134 -0.574844 -0.108574 0.440743 -0.112
#> marital_statusconsensual union -0.144878 -0.185376 -0.145006 -0.101527 -0.145
#> marital_statusNA 0.008810 -0.141012 0.009435 0.153243 0.017
#> study_years1 - 3 -0.247662 -0.319990 -0.247953 -0.169492 -0.245
#> study_years4 - 7 -0.299009 -0.368839 -0.298158 -0.229344 -0.298
#> study_years8 - 11 -0.248243 -0.322849 -0.248805 -0.175965 -0.249
#> study_years> 12 -0.091112 -0.189095 -0.091767 0.009038 -0.092
#> study_yearsNA -0.186889 -0.385197 -0.188717 0.018953 -0.171
#> consult_num1 - 3 -0.435789 -0.496850 -0.434688 -0.378305 -0.435
#> consult_num4 - 6 -0.642000 -0.704266 -0.642079 -0.582131 -0.639
#> consult_num> 7 -0.911590 -0.981541 -0.912101 -0.841647 -0.910
#> consult_numNA -0.041674 -0.192021 -0.040905 0.109398 -0.039
#> remoteness 0.183747 -0.089686 0.181960 0.461389 -0.106
#> prop_tap_toilet -0.256776 -0.821741 -0.266582 0.377560 -0.088
#> -
#> Acceptance probabilty:
#> Mean 2.5% 50% 97.5%
#> alpha 0.74479 0.03629 0.85692 1
#> -
#> Smooth terms:
#> Mean 2.5% 50% 97.5%
#> s(age).tau21 3.360e+00 8.782e-01 2.530e+00 1.048e+01
#> s(age).alpha 8.903e-01 2.618e-01 9.870e-01 1.000e+00
#> s(age).edf 5.681e+00 4.731e+00 5.654e+00 6.718e+00
#> s(wk_ini).tau21 3.868e+00 5.671e-01 2.882e+00 1.263e+01
#> s(wk_ini).alpha 9.698e-01 8.189e-01 9.942e-01 1.000e+00
#> s(wk_ini).edf 7.849e+00 6.338e+00 7.983e+00 8.687e+00
#> s(rivwk_conception).tau21 2.329e-04 3.466e-05 1.462e-04 9.944e-04
#> s(rivwk_conception).alpha 9.905e-01 9.398e-01 9.991e-01 1.000e+00
#> s(rivwk_conception).edf 3.640e+00 2.331e+00 3.525e+00 5.542e+00
#> s(res_muni).tau21 3.519e-02 2.148e-02 3.377e-02 5.598e-02
#> s(res_muni).alpha 7.801e-01 2.618e-01 8.676e-01 1.000e+00
#> s(res_muni).edf 3.796e+01 3.577e+01 3.802e+01 3.988e+01
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 3.151e-03 6.324e-05 9.774e-04 1.834e-02
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha 9.890e-01 9.252e-01 9.995e-01 1.000e+00
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 3.527e+00 2.060e+00 2.787e+00 8.722e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 4.219e-02 7.101e-04 3.435e-02 1.326e-01
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha 9.433e-01 6.942e-01 9.909e-01 1.000e+00
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 1.275e+01 2.861e+00 1.285e+01 2.221e+01
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 9.263e-03 7.401e-05 2.694e-03 5.700e-02
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha 9.540e-01 6.860e-01 9.974e-01 1.000e+00
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 5.598e+00 2.110e+00 4.497e+00 1.367e+01
#> parameters
#> s(age).tau21 0.668
#> s(age).alpha NA
#> s(age).edf 4.523
#> s(wk_ini).tau21 9.077
#> s(wk_ini).alpha NA
#> s(wk_ini).edf 8.586
#> s(rivwk_conception).tau21 0.000
#> s(rivwk_conception).alpha NA
#> s(rivwk_conception).edf 2.314
#> s(res_muni).tau21 8.008
#> s(res_muni).alpha NA
#> s(res_muni).edf 42.975
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).tau21 0.000
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).alpha NA
#> s(neg_mbsi_mean_1wk,pos_mbsi_mean_1wk).edf 1.756
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).tau21 0.045
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).alpha NA
#> s(neg_mckee_mean_8wk,pos_mckee_mean_8wk).edf 14.381
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).tau21 0.036
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).alpha NA
#> s(neg_ext_mbsi_mean_8wk,pos_ext_mbsi_mean_8wk).edf 12.031
#> ---
#> Sampler summary:
#> -
#> runtime = 24675.35
#> ---
#> Optimizer summary:
#> -
#> AICc = 120446.8 edf = 113.5654 logLik = -60109.78
#> logPost = -73466.1 nobs = 291479 runtime = 857.481Parametric 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
#> 7.083 0.250 7.408