Spatial 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: 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-17-sp-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 ~ s(longitud, latitud) + 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 5118459. logPost -3673047 logLik -2559209 edf 19.895 eps 0.2509 iteration 1"
#> [2] "AICc 4825471. logPost -2417561 logLik -2412722 edf 13.547 eps 0.1288 iteration 2"
#> [3] "AICc 4613570. logPost -2307717 logLik -2306732 edf 52.482 eps 0.0509 iteration 3"
#> [4] "AICc 4492049. logPost -2246453 logLik -2245952 edf 72.047 eps 0.0301 iteration 4"
#> [5] "AICc 4458524. logPost -2229670 logLik -2229189 edf 73.172 eps 0.0166 iteration 5"
#> [6] "AICc 4456442. logPost -2228628 logLik -2228147 edf 73.414 eps 0.0046 iteration 6"
#> [7] "AICc 4456434. logPost -2228624 logLik -2228143 edf 73.454 eps 0.0002 iteration 7"
#> [8] "AICc 4456434. logPost -2228624 logLik -2228143 edf 73.437 eps 0.0000 iteration 8"
#> [9] "AICc 4456434. logPost -2228624 logLik -2228143 edf 73.437 eps 0.0000 iteration 8"
#> [10] "elapsed time: 13.00sec"
#> [11] "Starting the sampler..."
#> [12] "Starting the sampler..."
#> [13] "Starting the sampler..."
#> [14] "Starting the sampler..."
#> [15] ""
#> [16] "| | 0% 11.31min"
#> [17] "| | 0% 11.44min"
#> [18] ""
#> [19] "|| | 0% | 0% 11.65min11.62min"
#> [20] "|* | 5% 10.63min 33.58sec"
#> [21] "|* | 5% 10.68min 33.71sec"
#> [22] ""
#> [23] "|* | 5% 10.96min 34.62sec|* | 5% 10.96min 34.60sec"
#> [24] "|** | 10% 9.79min 1.09min"
#> [25] "|** | 10% 9.81min 1.09min"
#> [26] "|** | 10% 10.02min 1.11min"
#> [27] "|** | 10% 10.31min 1.15min"
#> [28] "|*** | 15% 8.97min 1.58min"
#> [29] "|*** | 15% 8.98min 1.58min"
#> [30] "|*** | 15% 9.13min 1.61min"
#> [31] "|*** | 15% 9.59min 1.69min"
#> [32] "|**** | 20% 8.33min 2.08min"
#> [33] "|**** | 20% 8.35min 2.09min"
#> [34] "|**** | 20% 8.45min 2.11min"
#> [35] "|**** | 20% 8.77min 2.19min"
#> [36] "|***** | 25% 7.76min 2.59min"
#> [37] "|***** | 25% 7.77min 2.59min"
#> [38] "|***** | 25% 7.85min 2.62min"
#> [39] "|***** | 25% 8.10min 2.70min"
#> [40] "|****** | 30% 7.21min 3.09min"
#> [41] "|****** | 30% 7.22min 3.09min"
#> [42] "|****** | 30% 7.28min 3.12min"
#> [43] "|****** | 30% 7.47min 3.20min"
#> [44] "|******* | 35% 6.66min 3.59min"
#> [45] "|******* | 35% 6.67min 3.59min"
#> [46] "|******* | 35% 6.71min 3.61min"
#> [47] "|******* | 35% 6.87min 3.70min"
#> [48] "|******** | 40% 6.13min 4.08min"
#> [49] "|******** | 40% 6.14min 4.09min"
#> [50] "|******** | 40% 6.17min 4.11min"
#> [51] "|******** | 40% 6.30min 4.20min"
#> [52] "|********* | 45% 5.60min 4.59min"
#> [53] "|********* | 45% 5.62min 4.60min"
#> [54] "|********* | 45% 5.64min 4.61min"
#> [55] "|********* | 45% 5.74min 4.70min"
#> [56] "|********** | 50% 5.09min 5.09min"
#> [57] "|********** | 50% 5.09min 5.09min"
#> [58] "|********** | 50% 5.12min 5.12min"
#> [59] "|********** | 50% 5.19min 5.19min"
#> [60] "|*********** | 55% 4.57min 5.58min"
#> [61] "|*********** | 55% 4.58min 5.59min"
#> [62] "|*********** | 55% 4.59min 5.62min"
#> [63] "|*********** | 55% 4.66min 5.69min"
#> [64] "|************ | 60% 4.06min 6.09min"
#> [65] "|************ | 60% 4.06min 6.09min"
#> [66] "|************ | 60% 4.08min 6.11min"
#> [67] "|************ | 60% 4.13min 6.19min"
#> [68] "|************* | 65% 3.55min 6.59min"
#> [69] "|************* | 65% 3.55min 6.59min"
#> [70] "|************* | 65% 3.56min 6.61min"
#> [71] "|************* | 65% 3.61min 6.70min"
#> [72] "|************** | 70% 3.04min 7.08min"
#> [73] "|************** | 70% 3.04min 7.09min"
#> [74] "|************** | 70% 3.05min 7.11min"
#> [75] "|************** | 70% 3.08min 7.19min"
#> [76] "|*************** | 75% 2.53min 7.58min"
#> [77] "|*************** | 75% 2.53min 7.59min"
#> [78] "|*************** | 75% 2.54min 7.62min"
#> [79] "|*************** | 75% 2.56min 7.69min"
#> [80] "|**************** | 80% 2.02min 8.08min"
#> [81] "|**************** | 80% 2.02min 8.09min"
#> [82] "|**************** | 80% 2.03min 8.11min"
#> [83] "|**************** | 80% 2.05min 8.19min"
#> [84] "|***************** | 85% 1.51min 8.58min"
#> [85] "|***************** | 85% 1.51min 8.58min"
#> [86] "|***************** | 85% 1.52min 8.61min"
#> [87] "|***************** | 85% 1.53min 8.68min"
#> [88] "|****************** | 90% 1.01min 9.07min"
#> [89] "|****************** | 90% 1.01min 9.08min"
#> [90] "|****************** | 90% 1.01min 9.11min"
#> [91] "|****************** | 90% 1.02min 9.18min"
#> [92] "|******************* | 95% 30.25sec 9.58min"
#> [93] "|******************* | 95% 30.28sec 9.59min"
#> [94] "|******************* | 95% 30.34sec 9.61min"
#> [95] "|******************* | 95% 30.58sec 9.68min"
#> [96] "|********************| 100% 0.00sec 10.07min"
#> [97] ""
#> [98] "|********************| 100% 0.00sec 10.09min"
#> [99] ""
#> [100] "|********************| 100% 0.00sec 10.10min"
#> [101] ""
#> [102] "|********************| 100% 0.00sec 10.14min"system.time(saveRDS(bamlss_model, file = path_modelled_data))#> user system elapsed
#> 0.799 0.004 0.803saveRDS(form, file = path_modelled_form)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 2460.613 10.202 637.431