Full model with flexible prior: 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.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, "pre-12-full-re-p1.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 = 5, b = 1000)
form_aux <- premature ~ 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 168427.6 logPost -1626224 logLik -83963.4 edf 250.18 eps 0.5528 iteration 1"
#> [2] "AICc 155732.8 logPost -2630889 logLik -77588.0 edf 278.10 eps 0.2989 iteration 2"
#> [3] "AICc 153777.7 logPost -1162819 logLik -76602.2 edf 286.36 eps 0.1108 iteration 3"
#> [4] "AICc 153222.1 logPost -857626. logLik -76322.5 edf 288.27 eps 0.0550 iteration 4"
#> [5] "AICc 153055.5 logPost -755235. logLik -76239.0 edf 288.47 eps 0.0273 iteration 5"
#> [6] "AICc 153001.1 logPost -711519. logLik -76211.7 edf 288.54 eps 0.0137 iteration 6"
#> [7] "AICc 152978.8 logPost -694297. logLik -76200.6 edf 288.52 eps 0.0103 iteration 7"
#> [8] "AICc 152967.6 logPost -686469. logLik -76195.0 edf 288.49 eps 0.0046 iteration 8"
#> [9] "AICc 152961.1 logPost -682558. logLik -76191.8 edf 288.45 eps 0.0030 iteration 9"
#> [10] "AICc 152956.9 logPost -680459. logLik -76189.7 edf 288.42 eps 0.0020 iteration 10"
#> [11] "AICc 152954.0 logPost -679197. logLik -76188.3 edf 288.40 eps 0.0028 iteration 11"
#> [12] "AICc 152951.9 logPost -678320. logLik -76187.2 edf 288.39 eps 0.0014 iteration 12"
#> [13] "AICc 152950.3 logPost -677637. logLik -76186.5 edf 288.38 eps 0.0009 iteration 13"
#> [14] "AICc 152949.1 logPost -677068. logLik -76185.9 edf 288.37 eps 0.0008 iteration 14"
#> [15] "AICc 152948.2 logPost -676560. logLik -76185.4 edf 288.37 eps 0.0006 iteration 15"
#> [16] "AICc 152947.4 logPost -676096. logLik -76185.0 edf 288.37 eps 0.0005 iteration 16"
#> [17] "AICc 152946.8 logPost -675665. logLik -76184.7 edf 288.37 eps 0.0005 iteration 17"
#> [18] "AICc 152946.3 logPost -675260. logLik -76184.5 edf 288.37 eps 0.0006 iteration 18"
#> [19] "AICc 152945.9 logPost -674882. logLik -76184.3 edf 288.37 eps 0.0005 iteration 19"
#> [20] "AICc 152945.6 logPost -674525. logLik -76184.1 edf 288.37 eps 0.0006 iteration 20"
#> [21] "AICc 152945.3 logPost -674188. logLik -76184.0 edf 288.38 eps 0.0041 iteration 21"
#> [22] "AICc 152945.1 logPost -673872. logLik -76183.8 edf 288.38 eps 0.0003 iteration 22"
#> [23] "AICc 152944.9 logPost -673575. logLik -76183.7 edf 288.38 eps 0.0002 iteration 23"
#> [24] "AICc 152944.7 logPost -673299. logLik -76183.6 edf 288.39 eps 0.0002 iteration 24"
#> [25] "AICc 152944.5 logPost -673039. logLik -76183.6 edf 288.39 eps 0.0002 iteration 25"
#> [26] "AICc 152944.4 logPost -672800. logLik -76183.5 edf 288.40 eps 0.0001 iteration 26"
#> [27] "AICc 152944.3 logPost -672579. logLik -76183.4 edf 288.40 eps 0.0001 iteration 27"
#> [28] "AICc 152944.2 logPost -672375. logLik -76183.4 edf 288.40 eps 0.0001 iteration 28"
#> [29] "AICc 152944.1 logPost -672187. logLik -76183.3 edf 288.41 eps 0.0001 iteration 29"
#> [30] "AICc 152944.1 logPost -672014. logLik -76183.3 edf 288.41 eps 0.0001 iteration 30"
#> [31] "AICc 152944.0 logPost -671857. logLik -76183.3 edf 288.41 eps 0.0001 iteration 31"
#> [32] "AICc 152943.9 logPost -671709. logLik -76183.2 edf 288.42 eps 0.0001 iteration 32"
#> [33] "AICc 152943.9 logPost -671576. logLik -76183.2 edf 288.42 eps 0.0001 iteration 33"
#> [34] "AICc 152943.9 logPost -671457. logLik -76183.2 edf 288.42 eps 0.0001 iteration 34"
#> [35] "AICc 152943.8 logPost -671346. logLik -76183.2 edf 288.43 eps 0.0001 iteration 35"
#> [36] "AICc 152943.8 logPost -671247. logLik -76183.2 edf 288.43 eps 0.0000 iteration 36"
#> [37] "AICc 152943.8 logPost -671247. logLik -76183.2 edf 288.43 eps 0.0000 iteration 36"
#> [38] "elapsed time: 84.73min"
#> [39] "Starting the sampler..."
#> [40] "Starting the sampler..."
#> [41] "Starting the sampler..."
#> [42] "Starting the sampler..."
#> [43] ""
#> [44] "| | 0% 945.44min"
#> [45] "| | 0% 946.43min"
#> [46] "| | 0% 949.74min"
#> [47] "| | 0% 950.41min"
#> [48] "|* | 5% 882.67min 46.46min"
#> [49] "|* | 5% 882.83min 46.46min"
#> [50] "|* | 5% 886.32min 46.65min"
#> [51] "|* | 5% 887.51min 46.71min"
#> [52] "|** | 10% 835.96min 92.88min"
#> [53] "|** | 10% 836.61min 92.96min"
#> [54] "|** | 10% 839.36min 93.26min"
#> [55] "|** | 10% 839.79min 93.31min"
#> [56] "|*** | 15% 789.44min 139.31min"
#> [57] "|*** | 15% 791.28min 139.64min"
#> [58] "|*** | 15% 791.89min 139.74min"
#> [59] "|*** | 15% 792.59min 139.87min"
#> [60] "|**** | 20% 743.10min 185.77min"
#> [61] "|**** | 20% 744.74min 186.18min"
#> [62] "|**** | 20% 745.61min 186.40min"
#> [63] "|**** | 20% 745.65min 186.41min"
#> [64] "|***** | 25% 696.68min 232.23min"
#> [65] "|***** | 25% 698.43min 232.81min"
#> [66] "|***** | 25% 698.57min 232.86min"
#> [67] "|***** | 25% 699.24min 233.08min"
#> [68] "|****** | 30% 650.27min 278.69min"
#> [69] "|****** | 30% 652.01min 279.43min"
#> [70] "|****** | 30% 652.16min 279.50min"
#> [71] "|****** | 30% 652.29min 279.55min"
#> [72] "|******* | 35% 603.91min 325.18min"
#> [73] "|******* | 35% 605.43min 326.00min"
#> [74] "|******* | 35% 605.42min 326.00min"
#> [75] "|******* | 35% 605.71min 326.15min"
#> [76] "|******** | 40% 557.48min 371.65min"
#> [77] "|******** | 40% 558.74min 372.50min"
#> [78] "|******** | 40% 558.77min 372.51min"
#> [79] "|******** | 40% 559.20min 372.80min"
#> [80] "|********* | 45% 512.62min 419.41min"
#> [81] "|********* | 45% 513.97min 420.52min"
#> [82] "|********* | 45% 514.12min 420.65min"
#> [83] "|********* | 45% 514.51min 420.96min"
#> [84] "|********** | 50% 490.28min 490.28min"
#> [85] "|********** | 50% 491.48min 491.48min"
#> [86] "|********** | 50% 491.52min 491.52min"
#> [87] "|********** | 50% 491.80min 491.80min"
#> [88] "|*********** | 55% 449.60min 549.52min"
#> [89] "|*********** | 55% 450.31min 550.38min"
#> [90] "|*********** | 55% 450.42min 550.51min"
#> [91] "|*********** | 55% 450.52min 550.63min"
#> [92] "|************ | 60% 397.77min 596.66min"
#> [93] "|************ | 60% 398.43min 597.65min"
#> [94] "|************ | 60% 398.50min 597.75min"
#> [95] "|************ | 60% 398.59min 597.89min"
#> [96] "|************* | 65% 346.53min 643.55min"
#> [97] "|************* | 65% 347.04min 644.51min"
#> [98] "|************* | 65% 347.10min 644.61min"
#> [99] "|************* | 65% 347.14min 644.68min"
#> [100] "|************** | 70% 295.78min 690.16min"
#> [101] "|************** | 70% 296.16min 691.04min"
#> [102] "|************** | 70% 296.25min 691.25min"
#> [103] "|************** | 70% 296.25min 691.25min"
#> [104] "|*************** | 75% 245.52min 736.55min"
#> [105] "|*************** | 75% 245.90min 737.70min"
#> [106] "|*************** | 75% 245.94min 737.83min"
#> [107] "|*************** | 75% 245.96min 737.88min"
#> [108] "|**************** | 80% 195.77min 783.08min"
#> [109] "|**************** | 80% 196.07min 784.27min"
#> [110] "|**************** | 80% 196.11min 784.44min"
#> [111] "|**************** | 80% 196.11min 784.46min"
#> [112] "|***************** | 85% 146.42min 829.73min"
#> [113] "|***************** | 85% 146.62min 830.83min"
#> [114] "|***************** | 85% 146.65min 830.99min"
#> [115] "|***************** | 85% 146.66min 831.06min"
#> [116] "|****************** | 90% 97.55min 877.94min"
#> [117] "|****************** | 90% 97.66min 878.95min"
#> [118] "|****************** | 90% 97.68min 879.15min"
#> [119] "|****************** | 90% 97.69min 879.25min"
#> [120] "|******************* | 95% 48.71min 925.46min"
#> [121] "|******************* | 95% 48.76min 926.46min"
#> [122] "|******************* | 95% 48.77min 926.62min"
#> [123] "|******************* | 95% 48.78min 926.77min"
#> [124] "|********************| 100% 0.00sec 972.15min"
#> [125] ""
#> [126] "|********************| 100% 0.00sec 973.17min"
#> [127] ""
#> [128] "|********************| 100% 0.00sec 973.34min"
#> [129] ""
#> [130] "|********************| 100% 0.00sec 973.44min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 32.750 0.124 32.878saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 238907.411 136.272 63932.594