Growth model without rainfall: 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-25-growth-no-rain.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 <- 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")Run the model of interest and save results
{
sink(path_modelled_sink)
bamlss_model <- bamlss(
form, family = "binomial", data = bwdata_model,
n.iter = 8000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 132798.7 logPost -66661.2 logLik -66321.1 edf 78.229 eps 0.9950 iteration 1"
#> [2] "AICc 120871.3 logPost -60692.8 logLik -60353.0 edf 82.638 eps 0.2762 iteration 2"
#> [3] "AICc 120483.5 logPost -60481.1 logLik -60156.3 edf 85.368 eps 0.0524 iteration 3"
#> [4] "AICc 120471.5 logPost -60437.9 logLik -60150.3 edf 85.395 eps 0.0078 iteration 4"
#> [5] "AICc 120471.2 logPost -60436.2 logLik -60150.2 edf 85.383 eps 0.0017 iteration 5"
#> [6] "AICc 120471.2 logPost -60436.2 logLik -60150.2 edf 85.380 eps 0.0001 iteration 6"
#> [7] "AICc 120471.2 logPost -60436.2 logLik -60150.2 edf 85.379 eps 0.0000 iteration 7"
#> [8] "AICc 120471.2 logPost -60436.2 logLik -60150.2 edf 85.379 eps 0.0000 iteration 7"
#> [9] "elapsed time: 35.22sec"
#> [10] "Starting the sampler..."
#> [11] "Starting the sampler...Starting the sampler...Starting the sampler..."
#> [12] ""
#> [13] ""
#> [14] ""
#> [15] "| | 0% 175.47min"
#> [16] "| | 0% 178.92min"
#> [17] "| | 0% 179.00min"
#> [18] "| | 0% 179.42min"
#> [19] "|* | 5% 165.67min 8.72min"
#> [20] "|* | 5% 165.73min 8.72min"
#> [21] "|* | 5% 165.96min 8.73min"
#> [22] "|* | 5% 166.12min 8.74min"
#> [23] "|** | 10% 157.10min 17.46min"
#> [24] "|** | 10% 157.12min 17.46min"
#> [25] "|** | 10% 157.24min 17.47min"
#> [26] "|** | 10% 157.32min 17.48min"
#> [27] "|*** | 15% 148.57min 26.22min"
#> [28] "|*** | 15% 148.66min 26.23min"
#> [29] "|*** | 15% 148.71min 26.24min"
#> [30] "|*** | 15% 148.81min 26.26min"
#> [31] "|**** | 20% 139.67min 34.92min"
#> [32] "|**** | 20% 139.74min 34.93min"
#> [33] "|**** | 20% 139.76min 34.94min"
#> [34] "|**** | 20% 139.83min 34.96min"
#> [35] "|***** | 25% 130.92min 43.64min"
#> [36] ""
#> [37] "|***** | 25% 130.97min|***** | 25% 43.66min130.97min 43.66min"
#> [38] "|***** | 25% 131.03min 43.68min"
#> [39] "|****** | 30% 122.07min 52.32min"
#> [40] "|****** | 30% 122.09min 52.32min"
#> [41] "|****** | 30% 122.11min 52.33min"
#> [42] "|****** | 30% 122.16min 52.35min"
#> [43] "|******* | 35% 113.44min 61.08min"
#> [44] "|******* | 35% 113.47min 61.10min"
#> [45] "|******* | 35% 113.50min 61.11min"
#> [46] "|******* | 35% 113.54min 61.14min"
#> [47] "|******** | 40% 104.89min 69.93min"
#> [48] "|******** | 40% 104.90min 69.93min"
#> [49] "|******** | 40% 104.92min 69.95min"
#> [50] "|******** | 40% 104.94min 69.96min"
#> [51] "|********* | 45% 96.17min 78.69min"
#> [52] "|********* | 45% 96.19min 78.70min"
#> [53] "|********* | 45% 96.20min 78.71min"
#> [54] "|********* | 45% 96.22min 78.73min"
#> [55] "|********** | 50% 87.42min 87.42min"
#> [56] "|********** | 50% 87.43min 87.43min"
#> [57] "|********** | 50% 87.47min 87.47min"
#> [58] "|********** | 50% 87.47min 87.47min"
#> [59] "|*********** | 55% 78.68min 96.17min"
#> [60] "|*********** | 55% 78.70min 96.19min"
#> [61] "|*********** | 55% 78.75min 96.25min"
#> [62] "|*********** | 55% 78.74min 96.24min"
#> [63] "|************ | 60% 69.98min 104.97min"
#> [64] "|************ | 60% 70.00min 105.00min"
#> [65] "|************ | 60% 70.04min 105.06min"
#> [66] "|************ | 60% 70.05min 105.08min"
#> [67] "|************* | 65% 61.21min 113.68min"
#> [68] "|************* | 65% 61.25min 113.74min"
#> [69] "|************* | 65% 61.27min 113.79min"
#> [70] "|************* | 65% 61.30min 113.84min"
#> [71] "|************** | 70% 52.45min 122.38min"
#> [72] "|************** | 70% 52.48min 122.45min"
#> [73] "|************** | 70% 52.50min 122.51min"
#> [74] "|************** | 70% 52.53min 122.56min"
#> [75] "|*************** | 75% 43.70min 131.10min"
#> [76] "|*************** | 75% 43.73min 131.19min"
#> [77] "|*************** | 75% 43.75min 131.25min"
#> [78] "|*************** | 75% 43.77min 131.31min"
#> [79] "|**************** | 80% 34.96min 139.82min"
#> [80] "|**************** | 80% 34.98min 139.91min"
#> [81] "|**************** | 80% 34.99min 139.97min"
#> [82] "|**************** | 80% 35.01min 140.04min"
#> [83] "|***************** | 85% 26.21min 148.55min"
#> [84] "|***************** | 85% 26.23min 148.63min"
#> [85] "|***************** | 85% 26.24min 148.69min"
#> [86] "|***************** | 85% 26.26min 148.78min"
#> [87] "|****************** | 90% 17.47min 157.26min"
#> [88] "|****************** | 90% 17.48min 157.35min"
#> [89] "|****************** | 90% 17.49min 157.41min"
#> [90] "|****************** | 90% 17.50min 157.53min"
#> [91] "|******************* | 95% 8.74min 166.01min"
#> [92] "|******************* | 95% 8.74min 166.08min"
#> [93] "|******************* | 95% 8.75min 166.23min"
#> [94] "|******************* | 95% 8.75min 166.34min"
#> [95] "|********************| 100% 0.00sec 174.74min"
#> [96] ""
#> [97] "|********************| 100% 0.00sec 174.82min"
#> [98] ""
#> [99] "|********************| 100% 0.00sec 174.93min"
#> [100] ""
#> [101] "|********************| 100% 0.00sec 175.01min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 4.634 0.012 4.646saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 42074.587 23.962 10575.054