Linear remoteness with random effects and binning: fitting
Birthweight model including covariates at municipality level with linear effects.
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.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, "bw-muni-13-remote-re-bin.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 ~ remoteness + s(res_muni, bs = "re")
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, binning = TRUE,
n.iter = 1000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 4619787. logPost -3432973 logLik -2309881 edf 12.089 eps 0.0927 iteration 1"
#> [2] "AICc 4485810. logPost -2296595 logLik -2242879 edf 25.939 eps 0.0582 iteration 2"
#> [3] "AICc 4456869. logPost -2235069 logLik -2228391 edf 43.417 eps 0.0269 iteration 3"
#> [4] "AICc 4456382. logPost -2229063 logLik -2228145 edf 45.776 eps 0.0031 iteration 4"
#> [5] "AICc 4456379. logPost -2228509 logLik -2228143 edf 45.978 eps 0.0001 iteration 5"
#> [6] "AICc 4456379. logPost -2228489 logLik -2228143 edf 45.993 eps 0.0000 iteration 6"
#> [7] "AICc 4456379. logPost -2228489 logLik -2228143 edf 45.993 eps 0.0000 iteration 6"
#> [8] "elapsed time: 6.07sec"
#> [9] "Starting the sampler..."
#> [10] "Starting the sampler..."
#> [11] "Starting the sampler..."
#> [12] "Starting the sampler..."
#> [13] ""
#> [14] "| | 0% 9.18min"
#> [15] "| | 0% 9.28min"
#> [16] "| | 0% 9.45min"
#> [17] "| | 0% 9.45min"
#> [18] "|* | 5% 8.81min 27.81sec"
#> [19] "|* | 5% 8.83min 27.88sec"
#> [20] "|* | 5% 8.93min 28.21sec"
#> [21] "|* | 5% 8.96min 28.30sec"
#> [22] "|** | 10% 8.15min 54.34sec"
#> [23] "|** | 10% 8.16min 54.43sec"
#> [24] "|** | 10% 8.23min 54.89sec"
#> [25] "|** | 10% 8.26min 55.05sec"
#> [26] "|*** | 15% 7.58min 1.34min"
#> [27] "|*** | 15% 7.63min 1.35min"
#> [28] "|*** | 15% 7.63min 1.35min"
#> [29] "|*** | 15% 7.68min 1.35min"
#> [30] "|**** | 20% 7.09min 1.77min"
#> [31] "|**** | 20% 7.14min 1.78min"
#> [32] "|**** | 20% 7.17min 1.79min"
#> [33] "|**** | 20% 7.17min 1.79min"
#> [34] "|***** | 25% 6.62min 2.21min"
#> [35] "|***** | 25%"
#> [36] " 6.68min 2.23min|***** | 25% 6.68min 2.23min"
#> [37] "|***** | 25% 6.68min 2.23min"
#> [38] "|****** | 30% 6.19min 2.65min"
#> [39] "|****** | 30% 6.22min 2.67min"
#> [40] "|****** | 30% 6.22min 2.67min"
#> [41] "|****** | 30% 6.22min 2.67min"
#> [42] "|******* | 35% 5.73min 3.09min"
#> [43] "|******* | 35% 5.76min 3.10min"
#> [44] "|******* | 35% 5.76min 3.10min"
#> [45] "|******* | 35% 5.76min 3.10min"
#> [46] "|******** | 40% 5.30min 3.53min"
#> [47] "|******** | 40% 5.32min 3.55min"
#> [48] "|******** | 40% 5.32min 3.55min"
#> [49] "|******** | 40% 5.33min 3.55min"
#> [50] "|********* | 45% 4.88min 4.00min"
#> [51] "|********* | 45% 4.90min 4.01min"
#> [52] "|********* | 45% 4.91min 4.02min"
#> [53] "|********* | 45% 4.91min 4.02min"
#> [54] "|********** | 50% 4.43min 4.43min"
#> [55] "|********** | 50% 4.46min 4.46min"
#> [56] "|********** | 50% 4.46min 4.46min"
#> [57] "|********** | 50% 4.47min 4.47min"
#> [58] "|*********** | 55% 3.99min 4.88min"
#> [59] "|*********** | 55% 4.01min 4.90min"
#> [60] "|*********** | 55% 4.02min 4.91min"
#> [61] "|*********** | 55% 4.02min 4.92min"
#> [62] "|************ | 60% 3.56min 5.34min"
#> [63] "|************ | 60% 3.57min 5.36min"
#> [64] "|************ | 60% 3.58min 5.36min"
#> [65] "|************ | 60% 3.58min 5.37min"
#> [66] "|************* | 65% 3.12min 5.79min"
#> [67] "|************* | 65% 3.13min 5.81min"
#> [68] "|************* | 65% 3.13min 5.81min"
#> [69] "|************* | 65% 3.13min 5.81min"
#> [70] "|************** | 70% 2.67min 6.24min"
#> [71] "|************** | 70% 2.68min 6.25min"
#> [72] "|************** | 70% 2.68min 6.25min"
#> [73] "|************** | 70% 2.68min 6.26min"
#> [74] "|*************** | 75% 2.23min 6.69min"
#> [75] "|*************** | 75% 2.23min 6.70min"
#> [76] "|*************** | 75% 2.24min 6.71min"
#> [77] "|*************** | 75% 2.24min 6.72min"
#> [78] "|**************** | 80% 1.79min 7.15min"
#> [79] "|**************** | 80% 1.79min 7.16min"
#> [80] "|**************** | 80% 1.79min 7.16min"
#> [81] "|**************** | 80% 1.79min 7.17min"
#> [82] "|***************** | 85% 1.34min 7.60min"
#> [83] "|***************** | 85% 1.34min 7.62min"
#> [84] "|***************** | 85% 1.34min 7.62min"
#> [85] "|***************** | 85% 1.35min 7.63min"
#> [86] "|****************** | 90% 53.70sec 8.05min"
#> [87] "|****************** | 90% 53.75sec 8.06min"
#> [88] "|****************** | 90% 53.81sec 8.07min"
#> [89] "|****************** | 90% 53.83sec 8.07min"
#> [90] "|******************* | 95% 26.85sec 8.50min"
#> [91] "|******************* | 95% 26.87sec 8.51min"
#> [92] "|******************* | 95% 26.90sec 8.52min"
#> [93] "|******************* | 95% 26.94sec 8.53min"
#> [94] "|********************| 100% 0.00sec 8.95min"
#> [95] ""
#> [96] "|********************| 100% 0.00sec 8.95min"
#> [97] ""
#> [98] "|********************| 100% 0.00sec 8.96min"
#> [99] ""
#> [100] "|********************| 100% 0.00sec 8.97min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 0.831 0.004 0.835saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 2190.848 2.253 560.543