Full model with smoother 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.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-13-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 <- 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 -85039.2 logLik -83963.4 edf 250.18 eps 0.5528 iteration   1"
#>   [2] "AICc 155732.8 logPost -79068.1 logLik -77588.0 edf 278.10 eps 0.2989 iteration   2"
#>   [3] "AICc 153777.7 logPost -78171.3 logLik -76602.2 edf 286.36 eps 0.1108 iteration   3"
#>   [4] "AICc 153222.1 logPost -78057.3 logLik -76322.5 edf 288.27 eps 0.0550 iteration   4"
#>   [5] "AICc 153055.5 logPost -77884.7 logLik -76239.0 edf 288.47 eps 0.0273 iteration   5"
#>   [6] "AICc 153001.1 logPost -77797.2 logLik -76211.7 edf 288.54 eps 0.0137 iteration   6"
#>   [7] "AICc 152978.8 logPost -77781.6 logLik -76200.6 edf 288.52 eps 0.0103 iteration   7"
#>   [8] "AICc 152967.6 logPost -77776.4 logLik -76195.0 edf 288.49 eps 0.0046 iteration   8"
#>   [9] "AICc 152961.1 logPost -77773.2 logLik -76191.8 edf 288.45 eps 0.0030 iteration   9"
#>  [10] "AICc 152956.9 logPost -77771.0 logLik -76189.7 edf 288.42 eps 0.0020 iteration  10"
#>  [11] "AICc 152954.0 logPost -77769.5 logLik -76188.3 edf 288.40 eps 0.0028 iteration  11"
#>  [12] "AICc 152951.9 logPost -77768.3 logLik -76187.2 edf 288.39 eps 0.0014 iteration  12"
#>  [13] "AICc 152950.3 logPost -77767.5 logLik -76186.5 edf 288.38 eps 0.0009 iteration  13"
#>  [14] "AICc 152949.1 logPost -77766.8 logLik -76185.9 edf 288.37 eps 0.0008 iteration  14"
#>  [15] "AICc 152948.2 logPost -77766.3 logLik -76185.4 edf 288.37 eps 0.0006 iteration  15"
#>  [16] "AICc 152947.4 logPost -77766.0 logLik -76185.0 edf 288.37 eps 0.0005 iteration  16"
#>  [17] "AICc 152946.8 logPost -77765.7 logLik -76184.7 edf 288.37 eps 0.0005 iteration  17"
#>  [18] "AICc 152946.3 logPost -77765.4 logLik -76184.5 edf 288.37 eps 0.0006 iteration  18"
#>  [19] "AICc 152945.9 logPost -77765.3 logLik -76184.3 edf 288.37 eps 0.0005 iteration  19"
#>  [20] "AICc 152945.6 logPost -77765.1 logLik -76184.1 edf 288.37 eps 0.0006 iteration  20"
#>  [21] "AICc 152945.3 logPost -77765.0 logLik -76184.0 edf 288.38 eps 0.0041 iteration  21"
#>  [22] "AICc 152945.1 logPost -77764.9 logLik -76183.8 edf 288.38 eps 0.0003 iteration  22"
#>  [23] "AICc 152944.9 logPost -77764.9 logLik -76183.7 edf 288.38 eps 0.0002 iteration  23"
#>  [24] "AICc 152944.7 logPost -77764.8 logLik -76183.6 edf 288.39 eps 0.0002 iteration  24"
#>  [25] "AICc 152944.5 logPost -77764.8 logLik -76183.6 edf 288.39 eps 0.0002 iteration  25"
#>  [26] "AICc 152944.4 logPost -77764.8 logLik -76183.5 edf 288.40 eps 0.0001 iteration  26"
#>  [27] "AICc 152944.3 logPost -77764.8 logLik -76183.4 edf 288.40 eps 0.0001 iteration  27"
#>  [28] "AICc 152944.2 logPost -77764.7 logLik -76183.4 edf 288.40 eps 0.0001 iteration  28"
#>  [29] "AICc 152944.1 logPost -77764.7 logLik -76183.3 edf 288.41 eps 0.0001 iteration  29"
#>  [30] "AICc 152944.1 logPost -77764.7 logLik -76183.3 edf 288.41 eps 0.0001 iteration  30"
#>  [31] "AICc 152944.0 logPost -77764.7 logLik -76183.3 edf 288.41 eps 0.0001 iteration  31"
#>  [32] "AICc 152943.9 logPost -77764.7 logLik -76183.2 edf 288.42 eps 0.0001 iteration  32"
#>  [33] "AICc 152943.9 logPost -77764.7 logLik -76183.2 edf 288.42 eps 0.0001 iteration  33"
#>  [34] "AICc 152943.9 logPost -77764.7 logLik -76183.2 edf 288.42 eps 0.0001 iteration  34"
#>  [35] "AICc 152943.8 logPost -77764.7 logLik -76183.2 edf 288.43 eps 0.0001 iteration  35"
#>  [36] "AICc 152943.8 logPost -77764.7 logLik -76183.2 edf 288.43 eps 0.0000 iteration  36"
#>  [37] "AICc 152943.8 logPost -77764.7 logLik -76183.2 edf 288.43 eps 0.0000 iteration  36"
#>  [38] "elapsed time: 84.64min"                                                            
#>  [39] "Starting the sampler..."                                                           
#>  [40] "Starting the sampler..."                                                           
#>  [41] "Starting the sampler..."                                                           
#>  [42] "Starting the sampler..."                                                           
#>  [43] ""                                                                                  
#>  [44] "|                    |   0% 1049.17min"                                            
#>  [45] "|                    |   0% 1055.21min"                                            
#>  [46] "|                    |   0% 1034.40min"                                            
#>  [47] "|                    |   0% 1053.08min"                                            
#>  [48] "|*                   |   5% 1009.46min 53.13min"                                   
#>  [49] "|*                   |   5% 1011.91min 53.26min"                                   
#>  [50] "|*                   |   5% 1012.32min 53.28min"                                   
#>  [51] "|*                   |   5% 1015.38min 53.44min"                                   
#>  [52] "|**                  |  10% 920.93min 102.33min"                                   
#>  [53] "|**                  |  10% 921.74min 102.42min"                                   
#>  [54] "|**                  |  10% 921.53min 102.39min"                                   
#>  [55] "|**                  |  10% 923.63min 102.63min"                                   
#>  [56] "|***                 |  15% 843.83min 148.91min"                                   
#>  [57] "|***                 |  15% 844.50min 149.03min"                                   
#>  [58] "|***                 |  15% 844.50min 149.03min"                                   
#>  [59] "|***                 |  15% 846.85min 149.44min"                                   
#>  [60] "|****                |  20% 781.81min 195.45min"                                   
#>  [61] "|****                |  20% 782.50min 195.62min"                                   
#>  [62] "|****                |  20% 783.02min 195.75min"                                   
#>  [63] "|****                |  20% 784.71min 196.18min"                                   
#>  [64] "|*****               |  25% 725.75min 241.92min"                                   
#>  [65] "|*****               |  25% 726.01min 242.00min"                                   
#>  [66] "|*****               |  25% 727.39min 242.46min"                                   
#>  [67] "|*****               |  25% 728.48min 242.83min"                                   
#>  [68] "|******              |  30% 673.07min 288.46min"                                   
#>  [69] "|******              |  30% 673.06min 288.46min"                                   
#>  [70] "|******              |  30% 674.77min 289.19min"                                   
#>  [71] "|******              |  30% 675.05min 289.31min"                                   
#>  [72] "|*******             |  35% 621.90min 334.87min"                                   
#>  [73] "|*******             |  35% 622.54min 335.21min"                                   
#>  [74] "|*******             |  35% 623.39min 335.67min"                                   
#>  [75] "|*******             |  35% 623.99min 335.99min"                                   
#>  [76] "|********            |  40% 572.19min 381.46min"                                   
#>  [77] "|********            |  40% 572.70min 381.80min"                                   
#>  [78] "|********            |  40% 573.19min 382.13min"                                   
#>  [79] "|********            |  40% 573.73min 382.49min"                                   
#>  [80] "|*********           |  45% 523.19min 428.06min"                                   
#>  [81] "|*********           |  45% 523.59min 428.39min"                                   
#>  [82] "|*********           |  45% 523.92min 428.66min"                                   
#>  [83] "|*********           |  45% 524.22min 428.91min"                                   
#>  [84] "|**********          |  50% 474.56min 474.56min"                                   
#>  [85] "|**********          |  50% 474.96min 474.96min"                                   
#>  [86] "|**********          |  50% 475.26min 475.26min"                                   
#>  [87] "|**********          |  50% 475.32min 475.32min"                                   
#>  [88] "|***********         |  55% 426.24min 520.96min"                                   
#>  [89] "|***********         |  55% 426.75min 521.59min"                                   
#>  [90] "|***********         |  55% 426.94min 521.81min"                                   
#>  [91] "|***********         |  55% 427.07min 521.97min"                                   
#>  [92] "|************        |  60% 378.23min 567.35min"                                   
#>  [93] "|************        |  60% 378.74min 568.11min"                                   
#>  [94] "|************        |  60% 378.90min 568.34min"                                   
#>  [95] "|************        |  60% 379.13min 568.69min"                                   
#>  [96] "|*************       |  65% 330.46min 613.70min"                                   
#>  [97] "|*************       |  65% 330.93min 614.58min"                                   
#>  [98] "|*************       |  65% 331.12min 614.94min"                                   
#>  [99] "|*************       |  65% 331.38min 615.41min"                                   
#> [100] "|**************      |  70% 282.88min 660.06min"                                   
#> [101] "|**************      |  70% 283.32min 661.07min"                                   
#> [102] "|**************      |  70% 283.52min 661.54min"                                   
#> [103] "|**************      |  70% 283.75min 662.08min"                                   
#> [104] "|***************     |  75% 235.52min 706.57min"                                   
#> [105] "|***************     |  75% 235.85min 707.55min"                                   
#> [106] "|***************     |  75% 236.04min 708.11min"                                   
#> [107] "|***************     |  75% 236.21min 708.62min"                                   
#> [108] "|****************    |  80% 188.23min 752.91min"                                   
#> [109] "|****************    |  80% 188.51min 754.03min"                                   
#> [110] "|****************    |  80% 188.70min 754.79min"                                   
#> [111] "|****************    |  80% 188.81min 755.24min"                                   
#> [112] "|*****************   |  85% 141.07min 799.38min"                                   
#> [113] "|*****************   |  85% 141.25min 800.43min"                                   
#> [114] "|*****************   |  85% 141.43min 801.42min"                                   
#> [115] "|*****************   |  85% 141.50min 801.82min"                                   
#> [116] "|******************  |  90% 93.98min 845.79min"                                    
#> [117] "|******************  |  90% 94.11min 847.02min"                                    
#> [118] "|******************  |  90% 94.22min 847.95min"                                    
#> [119] "|******************  |  90% 94.28min 848.52min"                                    
#> [120] "|******************* |  95% 46.96min 892.32min"                                    
#> [121] "|******************* |  95% 47.03min 893.52min"                                    
#> [122] "|******************* |  95% 47.07min 894.38min"                                    
#> [123] "|******************* |  95% 47.12min 895.22min"                                    
#> [124] "|********************| 100%  0.00sec 940.05min"                                    
#> [125] ""                                                                                  
#> [126] "|********************| 100%  0.00sec 941.44min"                                    
#> [127] ""                                                                                  
#> [128] "|********************| 100%  0.00sec 941.77min"                                    
#> [129] ""                                                                                  
#> [130] "|********************| 100%  0.00sec 942.88min"
system.time(saveRDS(bamlss_model, file = path_modelled_data))
#>    user  system elapsed 
#>  33.037   0.172  33.215
saveRDS(form, file = path_modelled_form)

Time to execute the task

Only useful when executed with Rscript.

proc.time()
#>      user    system   elapsed 
#> 231667.26    111.62  62113.56