MRF model 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: MASSlibrary(sf)#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1library(dplyr)#>
#> Attaching package: 'dplyr'#> The following object is masked from 'package:MASS':
#>
#> select#> The following object is masked from 'package:bamlss':
#>
#> n#> The following object is masked from 'package:nlme':
#>
#> collapse#> The following objects are masked from 'package:stats':
#>
#> filter, lag#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, unionpath_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-19-mrf-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"))
K <- readRDS(file.path(path_processed, "ama_10_penalty.rds"))Define formula for our model
Now we define the same models as in the previous study.
form_sigma <- sigma ~ 1
form_mu <- born_weight ~ s(res_muni, bs = "mrf", xt = list("penalty" = K)) +
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 = 5000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 5119162. logPost -3736446 logLik -2559558 edf 23.352 eps 0.2579 iteration 1"
#> [2] "AICc 4825823. logPost -2418027 logLik -2412896 edf 14.893 eps 0.1308 iteration 2"
#> [3] "AICc 4614017. logPost -2307931 logLik -2306943 edf 64.966 eps 0.0509 iteration 3"
#> [4] "AICc 4492249. logPost -2246531 logLik -2246041 edf 83.555 eps 0.0301 iteration 4"
#> [5] "AICc 4458569. logPost -2229668 logLik -2229199 edf 85.273 eps 0.0167 iteration 5"
#> [6] "AICc 4456467. logPost -2228616 logLik -2228147 edf 85.748 eps 0.0046 iteration 6"
#> [7] "AICc 4456459. logPost -2228612 logLik -2228143 edf 85.817 eps 0.0002 iteration 7"
#> [8] "AICc 4456459. logPost -2228612 logLik -2228143 edf 85.767 eps 0.0000 iteration 8"
#> [9] "AICc 4456459. logPost -2228612 logLik -2228143 edf 85.767 eps 0.0000 iteration 8"
#> [10] "elapsed time: 49.23sec"
#> [11] "Starting the sampler..."
#> [12] "Starting the sampler..."
#> [13] "Starting the sampler..."
#> [14] "Starting the sampler..."
#> [15] ""
#> [16] "| | 0% 103.90min"
#> [17] "| | 0% 105.61min"
#> [18] "| | 0% 103.77min"
#> [19] "| | 0% 105.58min"
#> [20] "|* | 5% 105.05min 5.53min"
#> [21] "|* | 5% 105.81min 5.57min"
#> [22] "|* | 5% 105.90min 5.57min"
#> [23] "|* | 5% 106.00min 5.58min"
#> [24] "|** | 10% 98.71min 10.97min"
#> [25] "|** | 10% 100.16min 11.13min"
#> [26] "|** | 10% 100.20min 11.13min"
#> [27] "|** | 10% 100.21min 11.13min"
#> [28] "|*** | 15% 92.36min 16.30min"
#> [29] "|*** | 15% 93.32min 16.47min"
#> [30] "|*** | 15% 93.36min 16.48min"
#> [31] "|*** | 15% 93.38min 16.48min"
#> [32] "|**** | 20% 86.57min 21.64min"
#> [33] "|**** | 20% 87.25min 21.81min"
#> [34] "|**** | 20% 87.38min 21.85min"
#> [35] "|**** | 20% 87.49min 21.87min"
#> [36] "|***** | 25% 81.49min 27.16min"
#> [37] "|***** | 25% 82.13min 27.38min"
#> [38] "|***** | 25% 82.30min 27.43min"
#> [39] "|***** | 25% 82.60min 27.53min"
#> [40] "|****** | 30% 75.96min 32.56min"
#> [41] "|****** | 30% 76.56min 32.81min"
#> [42] "|****** | 30% 76.67min 32.86min"
#> [43] "|****** | 30% 76.92min 32.97min"
#> [44] "|******* | 35% 70.36min 37.89min"
#> [45] "|******* | 35% 70.79min 38.12min"
#> [46] "|******* | 35% 70.97min 38.22min"
#> [47] "|******* | 35% 71.31min 38.40min"
#> [48] "|******** | 40% 65.32min 43.55min"
#> [49] "|******** | 40% 65.63min 43.75min"
#> [50] "|******** | 40% 65.77min 43.84min"
#> [51] "|******** | 40% 65.98min 43.99min"
#> [52] "|********* | 45% 60.07min 49.14min"
#> [53] "|********* | 45% 60.34min 49.37min"
#> [54] "|********* | 45% 60.57min 49.56min"
#> [55] "|********* | 45% 60.78min 49.73min"
#> [56] "|********** | 50% 55.09min 55.09min"
#> [57] "|********** | 50% 55.30min 55.30min"
#> [58] "|********** | 50% 55.63min 55.63min"
#> [59] "|********** | 50% 55.79min 55.79min"
#> [60] "|*********** | 55% 49.65min 60.68min"
#> [61] "|*********** | 55% 49.78min 60.84min"
#> [62] "|*********** | 55% 50.08min 61.21min"
#> [63] "|*********** | 55% 50.19min 61.34min"
#> [64] "|************ | 60% 43.89min 65.84min"
#> [65] "|************ | 60% 43.98min 65.98min"
#> [66] "|************ | 60% 44.24min 66.36min"
#> [67] "|************ | 60% 44.33min 66.50min"
#> [68] "|************* | 65% 38.43min 71.38min"
#> [69] "|************* | 65% 38.55min 71.60min"
#> [70] "|************* | 65% 38.75min 71.96min"
#> [71] "|************* | 65% 38.84min 72.13min"
#> [72] "|************** | 70% 32.81min 76.56min"
#> [73] "|************** | 70% 32.90min 76.76min"
#> [74] "|************** | 70% 33.07min 77.17min"
#> [75] "|************** | 70% 33.16min 77.38min"
#> [76] "|*************** | 75% 27.22min 81.66min"
#> [77] "|*************** | 75% 27.30min 81.89min"
#> [78] "|*************** | 75% 27.44min 82.31min"
#> [79] "|*************** | 75% 27.50min 82.50min"
#> [80] "|**************** | 80% 21.71min 86.85min"
#> [81] "|**************** | 80% 21.78min 87.10min"
#> [82] "|**************** | 80% 21.91min 87.62min"
#> [83] "|**************** | 80% 21.95min 87.82min"
#> [84] "|***************** | 85% 16.28min 92.25min"
#> [85] "|***************** | 85% 16.33min 92.55min"
#> [86] "|***************** | 85% 16.43min 93.12min"
#> [87] "|***************** | 85% 16.45min 93.24min"
#> [88] "|****************** | 90% 10.86min 97.75min"
#> [89] "|****************** | 90% 10.90min 98.06min"
#> [90] "|****************** | 90% 10.96min 98.68min"
#> [91] "|****************** | 90% 10.96min 98.68min"
#> [92] "|******************* | 95% 5.43min 103.23min"
#> [93] "|******************* | 95% 5.45min 103.52min"
#> [94] "|******************* | 95% 5.48min 104.04min"
#> [95] "|******************* | 95% 5.48min 104.08min"
#> [96] "|********************| 100% 0.00sec 108.37min"
#> [97] ""
#> [98] "|********************| 100% 0.00sec 108.64min"
#> [99] ""
#> [100] "|********************| 100% 0.00sec 109.18min"
#> [101] ""
#> [102] "|********************| 100% 0.00sec 109.20min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 7.603 0.071 7.678saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 26016.188 13.049 6670.762