Growth model: fitting
Low birthweight model including covariates at municipality level with linear 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.constructpath_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")
path_modelled_data <- file.path(path_modelled, "lbw-20-growth-re.rds")
path_modelled_sink <- gsub("\\.rds$", "\\.txt", path_modelled_data)
path_modelled_form <- gsub("(\\.rds)$", "-form\\1", path_modelled_data)
bwdata_model <- fst::read_fst(file.path(path_processed, "bwdata_41_model.fst"))Define formula for our model
form_aux <- 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")
form <- update(form_aux, . ~ . +
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))Run the model of interest and save results
{
sink(path_modelled_sink)
bamlss_model <- bamlss(
form, family = "binomial", data = bwdata_model,
n.iter = 7000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 132902.9 logPost -66291.4 logLik -66306.5 edf 144.88 eps 0.7526 iteration 1"
#> [2] "AICc 120893.2 logPost -60643.4 logLik -60300.3 edf 146.16 eps 0.2720 iteration 2"
#> [3] "AICc 120464.9 logPost -60033.0 logLik -60106.8 edf 125.55 eps 0.0533 iteration 3"
#> [4] "AICc 120447.8 logPost -59873.1 logLik -60107.3 edf 116.48 eps 0.0102 iteration 4"
#> [5] "AICc 120447.3 logPost -59812.3 logLik -60108.9 edf 114.66 eps 0.0022 iteration 5"
#> [6] "AICc 120447.3 logPost -60031.4 logLik -60109.3 edf 114.20 eps 0.0008 iteration 6"
#> [7] "AICc 120447.1 logPost -62376.2 logLik -60109.5 edf 114.01 eps 0.0003 iteration 7"
#> [8] "AICc 120446.9 logPost -69473.7 logLik -60109.5 edf 113.82 eps 0.0007 iteration 8"
#> [9] "AICc 120446.8 logPost -72928.8 logLik -60109.6 edf 113.77 eps 0.0005 iteration 9"
#> [10] "AICc 120446.8 logPost -73412.3 logLik -60109.5 edf 113.78 eps 0.0003 iteration 10"
#> [11] "AICc 120446.8 logPost -73462.4 logLik -60109.5 edf 113.77 eps 0.0002 iteration 11"
#> [12] "AICc 120446.7 logPost -73467.2 logLik -60109.6 edf 113.74 eps 0.0002 iteration 12"
#> [13] "AICc 120446.7 logPost -73467.4 logLik -60109.6 edf 113.72 eps 0.0001 iteration 13"
#> [14] "AICc 120446.7 logPost -73467.2 logLik -60109.6 edf 113.69 eps 0.0001 iteration 14"
#> [15] "AICc 120446.7 logPost -73466.9 logLik -60109.6 edf 113.66 eps 0.0001 iteration 15"
#> [16] "AICc 120446.7 logPost -73466.7 logLik -60109.6 edf 113.63 eps 0.0001 iteration 16"
#> [17] "AICc 120446.7 logPost -73466.5 logLik -60109.7 edf 113.61 eps 0.0001 iteration 17"
#> [18] "AICc 120446.7 logPost -73466.3 logLik -60109.7 edf 113.59 eps 0.0001 iteration 18"
#> [19] "AICc 120446.7 logPost -73466.1 logLik -60109.7 edf 113.57 eps 0.0002 iteration 19"
#> [20] "AICc 120446.7 logPost -73466.0 logLik -60109.7 edf 113.56 eps 0.0000 iteration 20"
#> [21] "AICc 120446.7 logPost -73466.0 logLik -60109.7 edf 113.56 eps 0.0000 iteration 20"
#> [22] "elapsed time: 14.29min"
#> [23] "Starting the sampler...Starting the sampler...Starting the sampler..."
#> [24] "Starting the sampler..."
#> [25] ""
#> [26] ""
#> [27] ""
#> [28] ""
#> [29] ""
#> [30] ""
#> [31] "| | 0% 398.85min|| | | 0% | 0% 398.85min| 0% 398.85min398.85min"
#> [32] "|* | 5% 380.10min 20.01min"
#> [33] "|* | 5% 380.21min 20.01min"
#> [34] "|* | 5% 380.33min 20.02min"
#> [35] "|* | 5% 380.34min 20.02min"
#> [36] "|** | 10% 363.71min 40.41min"
#> [37] "|** | 10% 363.72min 40.41min"
#> [38] "|** | 10% 363.78min 40.42min"
#> [39] "|** | 10% 363.78min 40.42min"
#> [40] "|*** | 15% 344.79min 60.85min"
#> [41] "|*** | 15% 344.83min 60.85min"
#> [42] "|*** | 15% 344.87min 60.86min"
#> [43] "|*** | 15% 344.90min 60.87min"
#> [44] "|**** | 20% 327.01min 81.75min"
#> [45] "|**** | 20% 327.04min 81.76min"
#> [46] "|**** | 20% 327.06min 81.77min"
#> [47] "|**** | 20% 327.09min 81.77min"
#> [48] "|***** | 25% 307.94min 102.65min"
#> [49] "|***** | 25% 307.95min 102.65min"
#> [50] "|***** | 25% 307.97min 102.66min"
#> [51] "|***** | 25% 307.99min 102.66min"
#> [52] ""
#> [53] "|****** | 30% 287.90min 123.38min|****** | 30% 287.90min 123.38min"
#> [54] "|****** | 30% 287.93min 123.40min"
#> [55] "|****** | 30% 287.96min 123.41min"
#> [56] ""
#> [57] "|******* | 35% 267.54min 144.06min|******* | 35% 267.54min 144.06min"
#> [58] "|******* | 35% 267.56min 144.07min"
#> [59] "|******* | 35% 267.58min 144.08min"
#> [60] ""
#> [61] "|******** | 40% 246.92min 164.62min|******** | 40% 246.92min 164.62min"
#> [62] "|******** | 40% 246.94min 164.63min"
#> [63] "|******** | 40% 246.96min 164.64min"
#> [64] "|********* | 45% 226.42min 185.25min"
#> [65] "|********* | 45% 226.44min 185.27min"
#> [66] "|********* | 45% 226.51min 185.33min"
#> [67] "|********* | 45% 226.52min 185.33min"
#> [68] "|********** | 50% 205.68min 205.68min"
#> [69] "|********** | 50% 205.70min 205.70min"
#> [70] "|********** | 50% 205.76min 205.76min"
#> [71] "|********** | 50% 205.76min 205.76min"
#> [72] "|*********** | 55% 185.06min 226.19min"
#> [73] "|*********** | 55% 185.07min 226.20min"
#> [74] "|*********** | 55% 185.12min 226.26min"
#> [75] "|*********** | 55% 185.13min 226.27min"
#> [76] "|************ | 60% 164.47min 246.71min"
#> [77] "|************ | 60% 164.48min 246.72min"
#> [78] "|************ | 60% 164.52min 246.78min"
#> [79] "|************ | 60% 164.52min 246.79min"
#> [80] "|************* | 65% 143.87min 267.20min"
#> [81] "|************* | 65% 143.88min 267.21min"
#> [82] "|************* | 65% 143.91min 267.27min"
#> [83] "|************* | 65% 143.92min 267.27min"
#> [84] "|************** | 70% 123.30min 287.71min"
#> [85] "|************** | 70% 123.31min 287.72min"
#> [86] "|************** | 70% 123.33min 287.78min"
#> [87] "|************** | 70% 123.34min 287.79min"
#> [88] "|*************** | 75% 102.74min 308.21min"
#> [89] "|*************** | 75% 102.74min 308.22min"
#> [90] "|*************** | 75% 102.76min 308.28min"
#> [91] "|*************** | 75% 102.76min 308.29min"
#> [92] "|**************** | 80% 82.17min 328.69min"
#> [93] "|**************** | 80% 82.18min 328.70min"
#> [94] "|**************** | 80% 82.19min 328.76min"
#> [95] "|**************** | 80% 82.19min 328.77min"
#> [96] "|***************** | 85% 61.61min 349.14min"
#> [97] "|***************** | 85% 61.62min 349.16min"
#> [98] "|***************** | 85% 61.63min 349.21min"
#> [99] "|***************** | 85% 61.63min 349.22min"
#> [100] "|****************** | 90% 41.07min 369.62min"
#> [101] "|****************** | 90% 41.07min 369.63min"
#> [102] "|****************** | 90% 41.08min 369.69min"
#> [103] "|****************** | 90% 41.08min 369.70min"
#> [104] "|******************* | 95% 20.53min 390.14min"
#> [105] "|******************* | 95% 20.53min 390.15min"
#> [106] "|******************* | 95% 20.54min 390.20min"
#> [107] "|******************* | 95% 20.54min 390.21min"
#> [108] "|********************| 100% 0.00sec 410.63min"
#> [109] ""
#> [110] "|********************| 100% 0.00sec 410.64min"
#> [111] ""
#> [112] "|********************| 100% 0.00sec 410.69min"
#> [113] ""
#> [114] "|********************| 100% 0.00sec 410.69min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 6.985 0.036 7.072saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 99604.999 32.401 25622.346