MRF model with random effects: 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-18-mrf-re.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,
n.iter = 1000, burnin = 0, cores = 4, combine = FALSE, light = TRUE
)
sink()
}
readLines(path_modelled_sink)#> [1] "AICc 4465702. logPost -2455231 logLik -2232841 edf 9.6768 eps 0.0324 iteration 1"
#> [2] "AICc 4457087. logPost -2232089 logLik -2228510 edf 33.635 eps 0.0150 iteration 2"
#> [3] "AICc 4456437. logPost -2228674 logLik -2228149 edf 69.132 eps 0.0030 iteration 3"
#> [4] "AICc 4456449. logPost -2228525 logLik -2228144 edf 80.279 eps 0.0002 iteration 4"
#> [5] "AICc 4456449. logPost -2228522 logLik -2228144 edf 79.856 eps 0.0000 iteration 5"
#> [6] "AICc 4456449. logPost -2228522 logLik -2228144 edf 79.856 eps 0.0000 iteration 5"
#> [7] "elapsed time: 2.17min"
#> [8] "Starting the sampler..."
#> [9] "Starting the sampler..."
#> [10] "Starting the sampler..."
#> [11] "Starting the sampler..."
#> [12] ""
#> [13] "| | 0% 43.96min"
#> [14] "| | 0% 40.93min"
#> [15] "| | 0% 41.33min"
#> [16] "| | 0% 50.55min"
#> [17] "|* | 5% 41.77min 2.20min"
#> [18] "|* | 5% 41.12min 2.16min"
#> [19] "|* | 5% 43.04min 2.27min"
#> [20] "|* | 5% 44.43min 2.34min"
#> [21] "|** | 10% 40.29min 4.48min"
#> [22] "|** | 10% 41.15min 4.57min"
#> [23] "|** | 10% 41.28min 4.59min"
#> [24] "|** | 10% 43.22min 4.80min"
#> [25] "|*** | 15% 38.40min 6.78min"
#> [26] "|*** | 15% 39.08min 6.90min"
#> [27] "|*** | 15% 39.31min 6.94min"
#> [28] "|*** | 15% 40.90min 7.22min"
#> [29] "|**** | 20% 35.88min 8.97min"
#> [30] "|**** | 20% 36.18min 9.05min"
#> [31] "|**** | 20% 36.37min 9.09min"
#> [32] "|**** | 20% 37.32min 9.33min"
#> [33] "|***** | 25% 33.01min 11.00min"
#> [34] "|***** | 25% 33.27min 11.09min"
#> [35] "|***** | 25% 33.39min 11.13min"
#> [36] "|***** | 25% 34.14min 11.38min"
#> [37] "|****** | 30% 30.49min 13.07min"
#> [38] "|****** | 30% 30.66min 13.14min"
#> [39] "|****** | 30% 30.79min 13.19min"
#> [40] "|****** | 30% 31.32min 13.42min"
#> [41] "|******* | 35% 27.93min 15.04min"
#> [42] "|******* | 35% 28.06min 15.11min"
#> [43] "|******* | 35% 28.14min 15.15min"
#> [44] "|******* | 35% 28.60min 15.40min"
#> [45] "|******** | 40% 25.52min 17.01min"
#> [46] "|******** | 40% 25.63min 17.08min"
#> [47] "|******** | 40% 25.70min 17.13min"
#> [48] "|******** | 40% 26.07min 17.38min"
#> [49] "|********* | 45% 23.19min 18.97min"
#> [50] "|********* | 45% 23.29min 19.05min"
#> [51] "|********* | 45% 23.34min 19.10min"
#> [52] "|********* | 45% 23.66min 19.36min"
#> [53] "|********** | 50% 20.94min 20.94min"
#> [54] "|********** | 50% 21.01min 21.01min"
#> [55] "|********** | 50% 21.08min 21.08min"
#> [56] "|********** | 50% 21.32min 21.32min"
#> [57] "|*********** | 55% 18.71min 22.87min"
#> [58] "|*********** | 55% 18.77min 22.94min"
#> [59] "|*********** | 55% 18.82min 23.01min"
#> [60] "|*********** | 55% 19.01min 23.23min"
#> [61] "|************ | 60% 16.51min 24.76min"
#> [62] "|************ | 60% 16.55min 24.82min"
#> [63] "|************ | 60% 16.61min 24.91min"
#> [64] "|************ | 60% 16.76min 25.15min"
#> [65] "|************* | 65% 14.35min 26.64min"
#> [66] "|************* | 65% 14.38min 26.70min"
#> [67] "|************* | 65% 14.43min 26.81min"
#> [68] "|************* | 65% 14.56min 27.04min"
#> [69] "|************** | 70% 12.32min 28.74min"
#> [70] "|************** | 70% 12.35min 28.81min"
#> [71] "|************** | 70% 12.39min 28.92min"
#> [72] "|************** | 70% 12.50min 29.17min"
#> [73] "|*************** | 75% 10.28min 30.85min"
#> [74] "|*************** | 75% 10.30min 30.91min"
#> [75] "|*************** | 75% 10.36min 31.07min"
#> [76] "|*************** | 75% 10.45min 31.34min"
#> [77] "|**************** | 80% 8.24min 32.96min"
#> [78] "|**************** | 80% 8.26min 33.03min"
#> [79] "|**************** | 80% 8.30min 33.18min"
#> [80] "|**************** | 80% 8.36min 33.43min"
#> [81] "|***************** | 85% 6.20min 35.13min"
#> [82] "|***************** | 85% 6.21min 35.19min"
#> [83] "|***************** | 85% 6.24min 35.36min"
#> [84] "|***************** | 85% 6.28min 35.58min"
#> [85] "|****************** | 90% 4.12min 37.12min"
#> [86] "|****************** | 90% 4.13min 37.18min"
#> [87] "|****************** | 90% 4.15min 37.33min"
#> [88] "|****************** | 90% 4.17min 37.56min"
#> [89] "|******************* | 95% 2.06min 39.17min"
#> [90] "|******************* | 95% 2.07min 39.24min"
#> [91] "|******************* | 95% 2.07min 39.41min"
#> [92] "|******************* | 95% 2.09min 39.65min"
#> [93] "|********************| 100% 0.00sec 41.27min"
#> [94] ""
#> [95] "|********************| 100% 0.00sec 41.36min"
#> [96] ""
#> [97] "|********************| 100% 0.00sec 41.49min"
#> [98] ""
#> [99] "|********************| 100% 0.00sec 41.74min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 11.047 0.012 11.205saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 10249.420 34.947 2763.662