Full model with smoother prior: 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.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, "lbw-12-full-re-p2.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
hyper <- list(a = 10, b = 0.0001)
form_aux <- lbw ~ sex + born_race + birth_place +
marital_status + study_years + consult_num + s(age, xt = hyper) +
s(wk_ini, xt = hyper) + s(rivwk_conception, bs = "cc", xt = hyper) +
remoteness + prop_tap_toilet + s(res_muni, bs = "re", xt = hyper)
form <- update(form_aux, . ~ . +
s(neg_mbsi_mean_1wk, pos_mbsi_mean_1wk, k = 70, xt = hyper) +
s(neg_mckee_mean_8wk, pos_mckee_mean_8wk, k = 70, xt = hyper) +
s(neg_ext_mbsi_mean_8wk, pos_ext_mbsi_mean_8wk, k = 70, xt = hyper))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 139298.7 logPost -69820.0 logLik -69507.4 edf 141.88 eps 0.5483 iteration 1"
#> [2] "AICc 129929.2 logPost -66265.2 logLik -64781.3 edf 183.17 eps 0.2223 iteration 2"
#> [3] "AICc 129565.2 logPost -65649.0 logLik -64572.1 edf 210.29 eps 0.0380 iteration 3"
#> [4] "AICc 129537.6 logPost -65589.6 logLik -64550.5 edf 218.13 eps 0.0079 iteration 4"
#> [5] "AICc 129531.4 logPost -65618.6 logLik -64545.1 edf 220.39 eps 0.0029 iteration 5"
#> [6] "AICc 129529.0 logPost -65649.7 logLik -64543.1 edf 221.17 eps 0.0014 iteration 6"
#> [7] "AICc 129527.8 logPost -65662.1 logLik -64542.2 edf 221.48 eps 0.0008 iteration 7"
#> [8] "AICc 129527.2 logPost -65660.7 logLik -64541.8 edf 221.62 eps 0.0005 iteration 8"
#> [9] "AICc 129526.8 logPost -65654.4 logLik -64541.5 edf 221.68 eps 0.0003 iteration 9"
#> [10] "AICc 129526.5 logPost -65647.6 logLik -64541.3 edf 221.71 eps 0.0003 iteration 10"
#> [11] "AICc 129526.3 logPost -65641.9 logLik -64541.2 edf 221.73 eps 0.0002 iteration 11"
#> [12] "AICc 129526.2 logPost -65637.4 logLik -64541.2 edf 221.74 eps 0.0001 iteration 12"
#> [13] "AICc 129526.2 logPost -65634.1 logLik -64541.1 edf 221.75 eps 0.0001 iteration 13"
#> [14] "AICc 129526.1 logPost -65631.0 logLik -64541.1 edf 221.75 eps 0.0001 iteration 14"
#> [15] "AICc 129526.1 logPost -65629.0 logLik -64541.1 edf 221.76 eps 0.0001 iteration 15"
#> [16] "AICc 129526.1 logPost -65629.0 logLik -64541.1 edf 221.76 eps 0.0001 iteration 15"
#> [17] "elapsed time: 10.05min"
#> [18] "Starting the sampler..."
#> [19] "Starting the sampler..."
#> [20] "Starting the sampler..."
#> [21] "Starting the sampler..."
#> [22] ""
#> [23] "| | 0% 448.04min"
#> [24] "| | 0% 446.81min"
#> [25] "| | 0% 450.60min"
#> [26] "| | 0% 450.54min"
#> [27] "|* | 5% 419.67min 22.09min"
#> [28] ""
#> [29] "|* | 5% |* | 5%420.42min 420.36min22.13min 22.12min"
#> [30] "|* | 5% 420.34min 22.12min"
#> [31] "|** | 10% 397.61min 44.18min"
#> [32] ""
#> [33] "|** | 10% 398.08min 44.23min|** | 10% 398.11min 44.23min"
#> [34] "|** | 10% 398.07min 44.23min"
#> [35] "|*** | 15% 376.33min 66.41min"
#> [36] "|*** | 15% 376.70min 66.48min"
#> [37] "|*** | 15% 376.94min 66.52min"
#> [38] "|*** | 15% 376.96min 66.52min"
#> [39] "|**** | 20% 355.06min 88.76min"
#> [40] "|**** | 20% 355.26min 88.82min"
#> [41] "|**** | 20% 355.49min 88.87min"
#> [42] "|**** | 20% 355.47min 88.87min"
#> [43] "|***** | 25% 332.93min 110.98min"
#> [44] "|***** | 25% 333.08min 111.03min"
#> [45] "|***** | 25% 333.25min 111.08min"
#> [46] "|***** | 25% 333.24min 111.08min"
#> [47] "|****** | 30% 311.73min 133.60min"
#> [48] "|****** | 30% 311.99min 133.71min"
#> [49] "|****** | 30% 312.02min 133.72min"
#> [50] "|****** | 30% 312.16min 133.78min"
#> [51] "|******* | 35% 290.49min 156.42min"
#> [52] "|******* | 35% 290.69min 156.53min"
#> [53] "|******* | 35% 290.72min 156.54min"
#> [54] "|******* | 35% 290.83min 156.60min"
#> [55] "|******** | 40% 268.53min 179.02min"
#> [56] "|******** | 40% 268.70min 179.13min"
#> [57] "|******** | 40% 268.72min 179.15min"
#> [58] "|******** | 40% 268.81min 179.21min"
#> [59] "|********* | 45% 247.14min 202.21min"
#> [60] "|********* | 45% 247.23min 202.28min"
#> [61] "|********* | 45% 247.29min 202.33min"
#> [62] "|********* | 45% 247.32min 202.35min"
#> [63] "|********** | 50% 225.80min 225.80min"
#> [64] "|********** | 50% 225.87min 225.87min"
#> [65] "|********** | 50% 225.92min 225.92min"
#> [66] "|********** | 50% 225.96min 225.96min"
#> [67] "|*********** | 55% 204.05min 249.40min"
#> [68] "|*********** | 55% 204.11min 249.47min"
#> [69] "|*********** | 55% 204.15min 249.52min"
#> [70] "|*********** | 55% 204.18min 249.56min"
#> [71] "|************ | 60% 182.00min 272.99min"
#> [72] "|************ | 60% 182.04min 273.06min"
#> [73] "|************ | 60% 182.08min 273.11min"
#> [74] "|************ | 60% 182.10min 273.15min"
#> [75] "|************* | 65% 159.70min 296.58min"
#> [76] "|************* | 65% 159.73min 296.65min"
#> [77] "|************* | 65% 159.76min 296.70min"
#> [78] "|************* | 65% 159.79min 296.74min"
#> [79] "|************** | 70% 137.12min 319.95min"
#> [80] "|************** | 70% 137.15min 320.02min"
#> [81] "|************** | 70% 137.17min 320.06min"
#> [82] "|************** | 70% 137.19min 320.12min"
#> [83] "|*************** | 75% 114.31min 342.92min"
#> [84] "|*************** | 75% 114.33min 343.00min"
#> [85] "|*************** | 75% 114.35min 343.04min"
#> [86] "|*************** | 75% 114.37min 343.10min"
#> [87] "|**************** | 80% 91.49min 365.95min"
#> [88] "|**************** | 80% 91.51min 366.02min"
#> [89] "|**************** | 80% 91.52min 366.06min"
#> [90] "|**************** | 80% 91.53min 366.12min"
#> [91] "|***************** | 85% 68.63min 388.91min"
#> [92] "|***************** | 85% 68.64min 388.98min"
#> [93] "|***************** | 85% 68.65min 389.02min"
#> [94] "|***************** | 85% 68.66min 389.07min"
#> [95] "|****************** | 90% 45.76min 411.88min"
#> [96] "|****************** | 90% 45.77min 411.96min"
#> [97] "|****************** | 90% 45.78min 411.99min"
#> [98] "|****************** | 90% 45.78min 412.05min"
#> [99] "|******************* | 95% 22.89min 434.89min"
#> [100] "|******************* | 95% 22.89min 434.98min"
#> [101] "|******************* | 95% 22.89min 435.00min"
#> [102] "|******************* | 95% 22.90min 435.05min"
#> [103] "|********************| 100% 0.00sec 458.06min"
#> [104] ""
#> [105] "|********************| 100% 0.00sec 458.16min"
#> [106] ""
#> [107] "|********************| 100% 0.00sec 458.16min"
#> [108] ""
#> [109] "|********************| 100% 0.00sec 458.20min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 7.449 0.044 7.585saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 110725.012 43.991 28217.547