Gaussian full model: fitting
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.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")
path_modelled_data <- file.path(path_modelled, "bw-12-full-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_mu_aux <- born_weight ~ 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")
form_mu <- update(form_mu_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))
form_sigma <- sigma ~ sex + born_race + birth_place +
study_years + consult_num + s(age) +
s(wk_ini) + s(res_muni, bs = "re")
form <- list(form_mu, form_sigma)Run the model of interest and save results
{
sink(path_modelled_sink)
bamlss_model <- bamlss(
form, data = bwdata_model,
n.iter = 7000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 4445450. logPost -2449590 logLik -2222613 edf 111.37 eps 0.0369 iteration 1"
#> [2] "AICc 4438028. logPost -2222870 logLik -2218892 edf 121.73 eps 0.0081 iteration 2"
#> [3] "AICc 4434830. logPost -2219148 logLik -2217258 edf 156.55 eps 0.0046 iteration 3"
#> [4] "AICc 4433883. logPost -2218410 logLik -2216735 edf 205.85 eps 0.0029 iteration 4"
#> [5] "AICc 4433682. logPost -2218394 logLik -2216573 edf 267.39 eps 0.0011 iteration 5"
#> [6] "AICc 4433585. logPost -2218401 logLik -2216498 edf 293.90 eps 0.0007 iteration 6"
#> [7] "AICc 4433570. logPost -2218409 logLik -2216486 edf 298.16 eps 0.0002 iteration 7"
#> [8] "AICc 4433566. logPost -2218408 logLik -2216484 edf 298.54 eps 0.0000 iteration 8"
#> [9] "AICc 4433566. logPost -2218408 logLik -2216484 edf 298.54 eps 0.0000 iteration 8"
#> [10] "elapsed time: 21.40min"
#> [11] "Starting the sampler..."
#> [12] "Starting the sampler..."
#> [13] "Starting the sampler..."
#> [14] "Starting the sampler..."
#> [15] ""
#> [16] "| | 0% 1084.49min"
#> [17] "| | 0% 1070.54min"
#> [18] "| | 0% 1075.18min"
#> [19] "| | 0% 1089.47min"
#> [20] "|* | 5% 1003.19min 52.80min"
#> [21] "|* | 5% 1004.73min 52.88min"
#> [22] "|* | 5% 1005.43min 52.92min"
#> [23] "|* | 5% 1006.43min 52.97min"
#> [24] "|** | 10% 947.83min 105.31min"
#> [25] "|** | 10% 948.01min 105.33min"
#> [26] "|** | 10% 947.86min 105.32min"
#> [27] "|** | 10% 949.62min 105.51min"
#> [28] "|*** | 15% 880.17min 155.32min"
#> [29] "|*** | 15% 880.14min 155.32min"
#> [30] "|*** | 15% 880.44min 155.37min"
#> [31] "|*** | 15% 881.64min 155.58min"
#> [32] "|**** | 20% 817.21min 204.30min"
#> [33] "|**** | 20% 817.65min 204.41min"
#> [34] "|**** | 20% 817.74min 204.43min"
#> [35] "|**** | 20% 818.74min 204.69min"
#> [36] "|***** | 25% 758.52min 252.84min"
#> [37] "|***** | 25% 759.81min 253.27min"
#> [38] "|***** | 25% 760.28min 253.43min"
#> [39] "|***** | 25% 760.31min 253.44min"
#> [40] "|****** | 30% 703.65min 301.56min"
#> [41] "|****** | 30% 704.09min 301.75min"
#> [42] "|****** | 30% 704.76min 302.04min"
#> [43] "|****** | 30% 705.00min 302.14min"
#> [44] "|******* | 35% 650.37min 350.20min"
#> [45] "|******* | 35% 650.66min 350.35min"
#> [46] "|******* | 35% 651.28min 350.69min"
#> [47] "|******* | 35% 651.48min 350.80min"
#> [48] "|******** | 40% 598.09min 398.73min"
#> [49] "|******** | 40% 598.55min 399.03min"
#> [50] "|******** | 40% 598.91min 399.27min"
#> [51] "|******** | 40% 599.02min 399.34min"
#> [52] "|********* | 45% 547.34min 447.83min"
#> [53] "|********* | 45% 547.64min 448.07min"
#> [54] "|********* | 45% 548.15min 448.49min"
#> [55] "|********* | 45% 548.34min 448.64min"
#> [56] "|********** | 50% 496.34min 496.34min"
#> [57] "|********** | 50% 496.68min 496.68min"
#> [58] "|********** | 50% 496.89min 496.89min"
#> [59] "|********** | 50% 497.38min 497.38min"
#> [60] "|*********** | 55% 445.81min 544.88min"
#> [61] "|*********** | 55% 446.14min 545.28min"
#> [62] "|*********** | 55% 446.14min 545.28min"
#> [63] "|*********** | 55% 446.78min 546.06min"
#> [64] "|************ | 60% 395.72min 593.57min"
#> [65] "|************ | 60% 395.79min 593.69min"
#> [66] "|************ | 60% 395.91min 593.87min"
#> [67] "|************ | 60% 396.42min 594.63min"
#> [68] "|************* | 65% 345.73min 642.07min"
#> [69] "|************* | 65% 345.84min 642.27min"
#> [70] "|************* | 65% 345.94min 642.47min"
#> [71] "|************* | 65% 346.33min 643.19min"
#> [72] "|************** | 70% 295.95min 690.55min"
#> [73] "|************** | 70% 296.04min 690.76min"
#> [74] "|************** | 70% 296.17min 691.06min"
#> [75] "|************** | 70% 296.52min 691.88min"
#> [76] "|*************** | 75% 246.34min 739.02min"
#> [77] "|*************** | 75% 246.45min 739.35min"
#> [78] "|*************** | 75% 246.53min 739.60min"
#> [79] "|*************** | 75% 246.84min 740.53min"
#> [80] "|**************** | 80% 196.87min 787.48min"
#> [81] "|**************** | 80% 196.99min 787.97min"
#> [82] "|**************** | 80% 197.06min 788.26min"
#> [83] "|**************** | 80% 197.24min 788.97min"
#> [84] "|***************** | 85% 147.52min 835.96min"
#> [85] "|***************** | 85% 147.63min 836.60min"
#> [86] "|***************** | 85% 147.69min 836.93min"
#> [87] "|***************** | 85% 147.78min 837.43min"
#> [88] "|****************** | 90% 98.28min 884.49min"
#> [89] "|****************** | 90% 98.36min 885.23min"
#> [90] "|****************** | 90% 98.40min 885.56min"
#> [91] "|****************** | 90% 98.45min 886.03min"
#> [92] "|******************* | 95% 50.35min 956.58min"
#> [93] "|******************* | 95% 50.35min 956.58min"
#> [94] "|******************* | 95% 50.38min 957.14min"
#> [95] "|******************* | 95% 50.44min 958.28min"
#> [96] "|********************| 100% 0.00sec 1015.58min"
#> [97] ""
#> [98] "|********************| 100% 0.00sec 1015.66min"
#> [99] ""
#> [100] "|********************| 100% 0.00sec 1015.70min"
#> [101] ""
#> [102] "|********************| 100% 0.00sec 1016.54min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 58.989 0.316 59.310saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 245775.199 144.033 62851.567