Full model with centered covariates: 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-cd.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_aux <- premature ~ sex + born_race + birth_place +
marital_status + study_years + consult_num + s(age) +
s(wk_ini) + s(rivwk_conception, bs = "cc") +
remoteness_cd + prop_tap_toilet_cd + s(res_muni, bs = "re")
form <- update(form_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))Run the model of interest and save results
{
sink(path_modelled_sink)
bamlss_model <- bamlss(
form, family = "binomial", data = bwdata_model,
n.iter = 3000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 168427.6 logPost -84499.9 logLik -83963.4 edf 250.18 eps 0.5528 iteration 1"
#> [2] "AICc 155732.8 logPost -78372.0 logLik -77588.0 edf 278.10 eps 0.2989 iteration 2"
#> [3] "AICc 153777.7 logPost -77394.6 logLik -76602.2 edf 286.36 eps 0.1108 iteration 3"
#> [4] "AICc 153222.1 logPost -77222.7 logLik -76322.5 edf 288.27 eps 0.0550 iteration 4"
#> [5] "AICc 153055.5 logPost -77077.2 logLik -76239.0 edf 288.47 eps 0.0273 iteration 5"
#> [6] "AICc 153001.1 logPost -77012.1 logLik -76211.7 edf 288.54 eps 0.0137 iteration 6"
#> [7] "AICc 152978.8 logPost -76999.8 logLik -76200.6 edf 288.52 eps 0.0104 iteration 7"
#> [8] "AICc 152967.6 logPost -76994.6 logLik -76195.0 edf 288.49 eps 0.0046 iteration 8"
#> [9] "AICc 152961.1 logPost -76991.3 logLik -76191.8 edf 288.45 eps 0.0030 iteration 9"
#> [10] "AICc 152956.9 logPost -76989.2 logLik -76189.7 edf 288.42 eps 0.0020 iteration 10"
#> [11] "AICc 152954.0 logPost -76987.6 logLik -76188.3 edf 288.40 eps 0.0028 iteration 11"
#> [12] "AICc 152951.9 logPost -76986.5 logLik -76187.2 edf 288.39 eps 0.0014 iteration 12"
#> [13] "AICc 152950.3 logPost -76985.7 logLik -76186.5 edf 288.38 eps 0.0009 iteration 13"
#> [14] "AICc 152949.1 logPost -76985.0 logLik -76185.9 edf 288.37 eps 0.0008 iteration 14"
#> [15] "AICc 152948.2 logPost -76984.5 logLik -76185.4 edf 288.37 eps 0.0006 iteration 15"
#> [16] "AICc 152947.4 logPost -76984.1 logLik -76185.0 edf 288.37 eps 0.0005 iteration 16"
#> [17] "AICc 152946.8 logPost -76983.8 logLik -76184.7 edf 288.37 eps 0.0005 iteration 17"
#> [18] "AICc 152946.3 logPost -76983.6 logLik -76184.5 edf 288.37 eps 0.0006 iteration 18"
#> [19] "AICc 152945.9 logPost -76983.4 logLik -76184.3 edf 288.37 eps 0.0005 iteration 19"
#> [20] "AICc 152945.6 logPost -76983.3 logLik -76184.1 edf 288.37 eps 0.0006 iteration 20"
#> [21] "AICc 152945.3 logPost -76983.1 logLik -76184.0 edf 288.38 eps 0.0040 iteration 21"
#> [22] "AICc 152945.1 logPost -76983.0 logLik -76183.8 edf 288.38 eps 0.0003 iteration 22"
#> [23] "AICc 152944.9 logPost -76983.0 logLik -76183.7 edf 288.38 eps 0.0002 iteration 23"
#> [24] "AICc 152944.7 logPost -76982.9 logLik -76183.6 edf 288.39 eps 0.0002 iteration 24"
#> [25] "AICc 152944.5 logPost -76982.8 logLik -76183.6 edf 288.39 eps 0.0002 iteration 25"
#> [26] "AICc 152944.4 logPost -76982.8 logLik -76183.5 edf 288.40 eps 0.0003 iteration 26"
#> [27] "AICc 152944.3 logPost -76982.8 logLik -76183.4 edf 288.40 eps 0.0003 iteration 27"
#> [28] "AICc 152944.2 logPost -76982.8 logLik -76183.4 edf 288.41 eps 0.0001 iteration 28"
#> [29] "AICc 152944.2 logPost -76982.8 logLik -76183.4 edf 288.41 eps 0.0001 iteration 29"
#> [30] "AICc 152944.1 logPost -76982.7 logLik -76183.3 edf 288.41 eps 0.0002 iteration 30"
#> [31] "AICc 152944.0 logPost -76982.7 logLik -76183.3 edf 288.41 eps 0.0003 iteration 31"
#> [32] "AICc 152944.0 logPost -76982.7 logLik -76183.3 edf 288.42 eps 0.0001 iteration 32"
#> [33] "AICc 152943.9 logPost -76982.7 logLik -76183.2 edf 288.42 eps 0.0001 iteration 33"
#> [34] "AICc 152943.9 logPost -76982.7 logLik -76183.2 edf 288.42 eps 0.0001 iteration 33"
#> [35] "elapsed time: 77.58min"
#> [36] "Starting the sampler..."
#> [37] "Starting the sampler..."
#> [38] "Starting the sampler..."
#> [39] "Starting the sampler..."
#> [40] ""
#> [41] "| | 0% 354.30min"
#> [42] "| | 0% 353.28min"
#> [43] "| | 0% 355.52min"
#> [44] "| | 0% 354.65min"
#> [45] "|* | 5% 331.33min 17.44min"
#> [46] "|* | 5% 332.08min 17.48min"
#> [47] "|* | 5% 332.39min 17.49min"
#> [48] "|* | 5% 332.61min 17.51min"
#> [49] "|** | 10% 313.52min 34.84min"
#> [50] "|** | 10% 314.32min 34.92min"
#> [51] "|** | 10% 314.82min 34.98min"
#> [52] "|** | 10% 314.70min 34.97min"
#> [53] "|*** | 15% 295.91min 52.22min"
#> [54] "|*** | 15% 296.61min 52.34min"
#> [55] "|*** | 15% 297.19min 52.44min"
#> [56] "|*** | 15% 297.16min 52.44min"
#> [57] "|**** | 20% 278.81min 69.70min"
#> [58] "|**** | 20% 279.57min 69.89min"
#> [59] "|**** | 20% 279.79min 69.95min"
#> [60] "|**** | 20% 280.12min 70.03min"
#> [61] "|***** | 25% 261.44min 87.15min"
#> [62] "|***** | 25% 261.91min 87.30min"
#> [63] "|***** | 25% 262.06min 87.35min"
#> [64] "|***** | 25% 262.50min 87.50min"
#> [65] "|****** | 30% 243.99min 104.57min"
#> [66] "|****** | 30% 244.24min 104.67min"
#> [67] "|****** | 30% 244.44min 104.76min"
#> [68] "|****** | 30% 244.89min 104.95min"
#> [69] "|******* | 35% 226.62min 122.03min"
#> [70] "|******* | 35% 226.70min 122.07min"
#> [71] "|******* | 35% 226.91min 122.18min"
#> [72] "|******* | 35% 227.37min 122.43min"
#> [73] "|******** | 40% 209.19min 139.46min"
#> [74] "|******** | 40% 209.22min 139.48min"
#> [75] "|******** | 40% 209.40min 139.60min"
#> [76] "|******** | 40% 209.86min 139.90min"
#> [77] "|********* | 45% 191.71min 156.85min"
#> [78] "|********* | 45% 191.81min 156.94min"
#> [79] "|********* | 45% 191.89min 157.00min"
#> [80] "|********* | 45% 192.35min 157.38min"
#> [81] "|********** | 50% 174.26min 174.26min"
#> [82] "|********** | 50% 174.37min 174.37min"
#> [83] "|********** | 50% 174.41min 174.41min"
#> [84] "|********** | 50% 174.82min 174.82min"
#> [85] "|*********** | 55% 156.83min 191.68min"
#> [86] "|*********** | 55% 156.88min 191.74min"
#> [87] "|*********** | 55% 157.00min 191.88min"
#> [88] "|*********** | 55% 157.31min 192.27min"
#> [89] "|************ | 60% 139.39min 209.08min"
#> [90] "|************ | 60% 139.41min 209.11min"
#> [91] "|************ | 60% 139.56min 209.35min"
#> [92] "|************ | 60% 139.81min 209.72min"
#> [93] "|************* | 65% 121.96min 226.49min"
#> [94] "|************* | 65% 121.96min 226.49min"
#> [95] "|************* | 65% 122.13min 226.81min"
#> [96] "|************* | 65% 122.31min 227.15min"
#> [97] "|************** | 70% 104.52min 243.87min"
#> [98] "|************** | 70% 104.53min 243.91min"
#> [99] "|************** | 70% 104.69min 244.27min"
#> [100] "|************** | 70% 104.83min 244.60min"
#> [101] "|*************** | 75% 87.10min 261.31min"
#> [102] "|*************** | 75% 87.10min 261.31min"
#> [103] "|*************** | 75% 87.25min 261.74min"
#> [104] "|*************** | 75% 87.33min 261.99min"
#> [105] "|**************** | 80% 69.69min 278.75min"
#> [106] "|**************** | 80% 69.69min 278.75min"
#> [107] "|**************** | 80% 69.80min 279.21min"
#> [108] "|**************** | 80% 69.85min 279.41min"
#> [109] "|***************** | 85% 52.27min 296.19min"
#> [110] "|***************** | 85% 52.27min 296.20min"
#> [111] "|***************** | 85% 52.34min 296.58min"
#> [112] "|***************** | 85% 52.39min 296.86min"
#> [113] "|****************** | 90% 34.85min 313.65min"
#> [114] "|****************** | 90% 34.85min 313.66min"
#> [115] "|****************** | 90% 34.89min 313.99min"
#> [116] "|****************** | 90% 34.92min 314.26min"
#> [117] "|******************* | 95% 17.43min 331.09min"
#> [118] "|******************* | 95% 17.43min 331.15min"
#> [119] "|******************* | 95% 17.44min 331.40min"
#> [120] "|******************* | 95% 17.46min 331.67min"
#> [121] "|********************| 100% 0.00sec 348.54min"
#> [122] ""
#> [123] "|********************| 100% 0.00sec 348.58min"
#> [124] ""
#> [125] "|********************| 100% 0.00sec 348.81min"
#> [126] ""
#> [127] "|********************| 100% 0.00sec 349.13min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 23.780 0.064 23.772saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 88968.454 59.667 26031.531