Full model without rainfall: fitting
Prematurity 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.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")
path_modelled_data <- file.path(path_modelled, "pre-15-full-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 <- 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")
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 171382.6 logPost -85885.4 logLik -85609.9 edf 81.365 eps 0.5494 iteration 1"
#> [2] "AICc 162371.6 logPost -81370.2 logLik -81100.5 edf 85.282 eps 0.1913 iteration 2"
#> [3] "AICc 161244.9 logPost -80856.5 logLik -80536.1 edf 86.241 eps 0.0598 iteration 3"
#> [4] "AICc 161037.1 logPost -80718.3 logLik -80432.4 edf 86.039 eps 0.0236 iteration 4"
#> [5] "AICc 160996.8 logPost -80690.3 logLik -80412.5 edf 85.835 eps 0.0103 iteration 5"
#> [6] "AICc 160988.9 logPost -80686.2 logLik -80408.7 edf 85.754 eps 0.0046 iteration 6"
#> [7] "AICc 160987.4 logPost -80685.4 logLik -80407.9 edf 85.720 eps 0.0020 iteration 7"
#> [8] "AICc 160987.1 logPost -80685.3 logLik -80407.8 edf 85.706 eps 0.0009 iteration 8"
#> [9] "AICc 160987.1 logPost -80685.3 logLik -80407.8 edf 85.699 eps 0.0004 iteration 9"
#> [10] "AICc 160987.1 logPost -80685.2 logLik -80407.8 edf 85.696 eps 0.0001 iteration 10"
#> [11] "AICc 160987.1 logPost -80685.2 logLik -80407.8 edf 85.695 eps 0.0000 iteration 11"
#> [12] "AICc 160987.1 logPost -80685.2 logLik -80407.8 edf 85.695 eps 0.0000 iteration 11"
#> [13] "elapsed time: 51.03sec"
#> [14] "Starting the sampler...Starting the sampler..."
#> [15] ""
#> [16] "Starting the sampler..."
#> [17] "Starting the sampler..."
#> [18] ""
#> [19] ""
#> [20] ""
#> [21] "|| | | 0% | 0% 169.80min169.72min | 0% 169.61min"
#> [22] "| | 0% 169.92min"
#> [23] "|* | 5% 159.55min 8.40min"
#> [24] ""
#> [25] "|* | 5% |* | 5%159.94min 159.94min 8.42min 8.42min"
#> [26] "|* | 5% 159.94min 8.42min"
#> [27] "|** | 10% 151.15min 16.79min"
#> [28] "|** | 10% 151.15min 16.79min"
#> [29] "|** | 10% 151.33min 16.81min"
#> [30] "|** | 10% 151.33min 16.81min"
#> [31] "|*** | 15% 142.85min 25.21min"
#> [32] "|*** | 15% 142.95min 25.23min"
#> [33] "|*** | 15% 143.04min 25.24min"
#> [34] "|*** | 15% 143.05min 25.24min"
#> [35] "|**** | 20% 134.41min 33.60min"
#> [36] "|**** | 20% 134.50min 33.63min"
#> [37] "|**** | 20% 134.57min 33.64min"
#> [38] "|**** | 20% 134.58min 33.65min"
#> [39] "|***** | 25% 126.06min 42.02min"
#> [40] "|***** | 25% 126.25min 42.08min"
#> [41] "|***** | 25% 126.29min 42.10min"
#> [42] "|***** | 25% 126.31min 42.10min"
#> [43] "|****** | 30% 117.78min 50.48min"
#> [44] "|****** | 30% 117.95min 50.55min"
#> [45] "|****** | 30% 117.98min 50.56min"
#> [46] "|****** | 30% 118.00min 50.57min"
#> [47] "|******* | 35% 109.40min 58.91min"
#> [48] "|******* | 35% 109.59min 59.01min"
#> [49] "|******* | 35% 109.59min 59.01min"
#> [50] "|******* | 35% 109.63min 59.03min"
#> [51] "|******** | 40% 100.99min 67.33min"
#> [52] ""
#> [53] "|******** | 40% 101.18min 67.45min|******** | 40% 101.18min 67.45min"
#> [54] "|******** | 40% 101.21min 67.47min"
#> [55] "|********* | 45% 92.55min 75.72min"
#> [56] "|********* | 45% 92.70min 75.84min"
#> [57] "|********* | 45% 92.71min 75.85min"
#> [58] "|********* | 45% 92.73min 75.87min"
#> [59] "|********** | 50% 84.24min 84.24min"
#> [60] "|********** | 50% 84.36min 84.36min"
#> [61] "|********** | 50% 84.37min 84.37min"
#> [62] "|********** | 50% 84.39min 84.39min"
#> [63] "|*********** | 55% 75.86min 92.72min"
#> [64] "|*********** | 55% 75.96min 92.84min"
#> [65] "|*********** | 55% 75.98min 92.86min"
#> [66] "|*********** | 55% 75.99min 92.88min"
#> [67] "|************ | 60% 67.47min 101.20min"
#> [68] "|************ | 60% 67.54min 101.31min"
#> [69] "|************ | 60% 67.56min 101.34min"
#> [70] "|************ | 60% 67.58min 101.37min"
#> [71] "|************* | 65% 59.05min 109.66min"
#> [72] "|************* | 65% 59.10min 109.76min"
#> [73] "|************* | 65% 59.13min 109.81min"
#> [74] "|************* | 65% 59.15min 109.85min"
#> [75] "|************** | 70% 50.63min 118.14min"
#> [76] "|************** | 70% 50.68min 118.25min"
#> [77] "|************** | 70% 50.70min 118.30min"
#> [78] "|************** | 70% 50.71min 118.33min"
#> [79] "|*************** | 75% 42.20min 126.61min"
#> [80] "|*************** | 75% 42.24min 126.73min"
#> [81] "|*************** | 75% 42.27min 126.81min"
#> [82] "|*************** | 75% 42.28min 126.83min"
#> [83] "|**************** | 80% 33.77min 135.08min"
#> [84] "|**************** | 80% 33.80min 135.21min"
#> [85] "|**************** | 80% 33.82min 135.30min"
#> [86] "|**************** | 80% 33.83min 135.31min"
#> [87] "|***************** | 85% 25.33min 143.55min"
#> [88] "|***************** | 85% 25.36min 143.69min"
#> [89] "|***************** | 85% 25.37min 143.78min"
#> [90] "|***************** | 85% 25.37min 143.78min"
#> [91] "|****************** | 90% 16.90min 152.06min"
#> [92] "|****************** | 90% 16.91min 152.18min"
#> [93] "|****************** | 90% 16.92min 152.28min"
#> [94] "|****************** | 90% 16.92min 152.29min"
#> [95] "|******************* | 95% 8.45min 160.53min"
#> [96] "|******************* | 95% 8.46min 160.65min"
#> [97] "|******************* | 95% 8.46min 160.76min"
#> [98] "|******************* | 95% 8.46min 160.78min"
#> [99] "|********************| 100% 0.00sec 168.96min"
#> [100] ""
#> [101] "|********************| 100% 0.00sec 169.05min"
#> [102] ""
#> [103] "|********************| 100% 0.00sec 169.15min"
#> [104] ""
#> [105] "|********************| 100% 0.00sec 169.17min"
system.time(saveRDS(bamlss_model, file = path_modelled_data))
#> user system elapsed
#> 4.570 0.016 4.586
saveRDS(form, file = path_modelled_form)
Time to execute the task
Only useful when executed with Rscript
.
proc.time()
#> user system elapsed
#> 40699.438 23.087 10240.537