Parametric remoteness: fitting


Birthweight model including parametric remoteness effect.

Load packages, read data and source custom scripts

Paths are defined relative to the git repository location.

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, "bw-muni-06-remote-param.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

Now we define the same models as in the previous study.

form_sigma <- sigma ~ 1

form_mu <- born_weight ~ poly(remoteness, 2)

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 = 1000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
    )
    sink()
}
readLines(path_modelled_sink)
#>  [1] "AICc 4458361. logPost -2229270 logLik -2229176 edf 4.0000 eps 0.0299 iteration   1"
#>  [2] "AICc 4458361. logPost -2229270 logLik -2229176 edf 4.0000 eps 0.0000 iteration   2"
#>  [3] "AICc 4458361. logPost -2229270 logLik -2229176 edf 4.0000 eps 0.0000 iteration   2"
#>  [4] "elapsed time:  0.44sec"                                                            
#>  [5] "Starting the sampler..."                                                           
#>  [6] "Starting the sampler..."                                                           
#>  [7] "Starting the sampler..."                                                           
#>  [8] "Starting the sampler..."                                                           
#>  [9] ""                                                                                  
#> [10] "|                    |   0%  5.69min"                                              
#> [11] "|                    |   0%  5.78min"                                              
#> [12] ""                                                                                  
#> [13] "|                    |   0%  5.70min|                    |   0%  5.70min"          
#> [14] "|*                   |   5%  4.42min 13.95sec"                                     
#> [15] "|*                   |   5%  4.41min 13.93sec"                                     
#> [16] "|*                   |   5%  4.43min 13.98sec"                                     
#> [17] "|*                   |   5%  4.57min 14.45sec"                                     
#> [18] "|**                  |  10%  4.07min 27.11sec"                                     
#> [19] "|**                  |  10%  4.07min 27.11sec"                                     
#> [20] "|**                  |  10%  4.11min 27.42sec"                                     
#> [21] "|**                  |  10%  4.15min 27.67sec"                                     
#> [22] "|***                 |  15%  3.81min 40.36sec"                                     
#> [23] "|***                 |  15%  3.81min 40.35sec"                                     
#> [24] "|***                 |  15%  3.81min 40.37sec"                                     
#> [25] "|***                 |  15%  3.87min 40.95sec"                                     
#> [26] "|****                |  20%  3.57min 53.52sec"                                     
#> [27] "|****                |  20%  3.57min 53.59sec"                                     
#> [28] "|****                |  20%  3.57min 53.60sec"                                     
#> [29] "|****                |  20%  3.61min 54.18sec"                                     
#> [30] "|*****               |  25%  3.33min  1.11min"                                     
#> [31] "|*****               |  25%  3.34min  1.11min"                                     
#> [32] "|*****               |  25%  3.34min  1.11min"                                     
#> [33] "|*****               |  25%  3.37min  1.12min"                                     
#> [34] "|******              |  30%  3.11min  1.33min"                                     
#> [35] "|******              |  30%  3.12min  1.34min"                                     
#> [36] "|******              |  30%  3.12min  1.34min"                                     
#> [37] "|******              |  30%  3.15min  1.35min"                                     
#> [38] "|*******             |  35%  2.89min  1.56min"                                     
#> [39] "|*******             |  35%  2.90min  1.56min"                                     
#> [40] "|*******             |  35%  2.91min  1.57min"                                     
#> [41] "|*******             |  35%  2.91min  1.57min"                                     
#> [42] "|********            |  40%  2.67min  1.78min"                                     
#> [43] "|********            |  40%  2.68min  1.79min"                                     
#> [44] "|********            |  40%  2.69min  1.79min"                                     
#> [45] "|********            |  40%  2.69min  1.79min"                                     
#> [46] "|*********           |  45%  2.46min  2.01min"                                     
#> [47] "|*********           |  45%  2.46min  2.01min"                                     
#> [48] "|*********           |  45%  2.46min  2.02min"                                     
#> [49] "|*********           |  45%  2.47min  2.02min"                                     
#> [50] "|**********          |  50%  2.24min  2.24min"                                     
#> [51] "|**********          |  50%  2.24min  2.24min"                                     
#> [52] "|**********          |  50%  2.25min  2.25min"                                     
#> [53] "|**********          |  50%  2.25min  2.25min"                                     
#> [54] "|***********         |  55%  2.01min  2.46min"                                     
#> [55] "|***********         |  55%  2.02min  2.47min"                                     
#> [56] "|***********         |  55%  2.02min  2.47min"                                     
#> [57] "|***********         |  55%  2.03min  2.48min"                                     
#> [58] "|************        |  60%  1.79min  2.69min"                                     
#> [59] "|************        |  60%  1.80min  2.69min"                                     
#> [60] "|************        |  60%  1.80min  2.70min"                                     
#> [61] "|************        |  60%  1.80min  2.71min"                                     
#> [62] "|*************       |  65%  1.57min  2.92min"                                     
#> [63] "|*************       |  65%  1.57min  2.92min"                                     
#> [64] "|*************       |  65%  1.57min  2.92min"                                     
#> [65] "|*************       |  65%  1.58min  2.93min"                                     
#> [66] "|**************      |  70%  1.35min  3.14min"                                     
#> [67] "|**************      |  70%  1.35min  3.14min"                                     
#> [68] "|**************      |  70%  1.35min  3.15min"                                     
#> [69] "|**************      |  70%  1.35min  3.16min"                                     
#> [70] "|***************     |  75%  1.12min  3.37min"                                     
#> [71] "|***************     |  75%  1.12min  3.37min"                                     
#> [72] "|***************     |  75%  1.12min  3.37min"                                     
#> [73] "|***************     |  75%  1.13min  3.39min"                                     
#> [74] "|****************    |  80% 53.96sec  3.60min"                                     
#> [75] "|****************    |  80% 53.99sec  3.60min"                                     
#> [76] "|****************    |  80% 53.99sec  3.60min"                                     
#> [77] "|****************    |  80% 54.22sec  3.61min"                                     
#> [78] "|*****************   |  85% 40.43sec  3.82min"                                     
#> [79] "|*****************   |  85% 40.46sec  3.82min"                                     
#> [80] "|*****************   |  85% 40.53sec  3.83min"                                     
#> [81] "|*****************   |  85% 40.66sec  3.84min"                                     
#> [82] "|******************  |  90% 26.95sec  4.04min"                                     
#> [83] "|******************  |  90% 26.98sec  4.05min"                                     
#> [84] "|******************  |  90% 27.02sec  4.05min"                                     
#> [85] "|******************  |  90% 27.11sec  4.07min"                                     
#> [86] "|******************* |  95% 13.49sec  4.27min"                                     
#> [87] "|******************* |  95% 13.51sec  4.28min"                                     
#> [88] "|******************* |  95% 13.52sec  4.28min"                                     
#> [89] "|******************* |  95% 13.57sec  4.30min"                                     
#> [90] "|********************| 100%  0.00sec  4.50min"                                     
#> [91] ""                                                                                  
#> [92] "|********************| 100%  0.00sec  4.50min"                                     
#> [93] ""                                                                                  
#> [94] "|********************| 100%  0.00sec  4.51min"                                     
#> [95] ""                                                                                  
#> [96] "|********************| 100%  0.00sec  4.51min"
system.time(saveRDS(bamlss_model, file = path_modelled_data))
#>    user  system elapsed 
#>   0.209   0.000   0.210
saveRDS(form, file = path_modelled_form)

Time to execute the task

Only useful when executed with Rscript.

proc.time()
#>     user   system  elapsed 
#> 1112.623    7.388  287.744